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Previous  studies  predicted  the  F-15B  high  angle  of 
attack  and  flat  spin  behavior  using  bifurcation  analysis. 
These  studies  varied  control  surface  deflections  to  find 
equilibrium  and  periodic  solutions.  The  purpose  of  this 

research  was  to  use  bifurcation  analysis  to  predict 
the  F-15B  high  angle  of  attack  and  flat  spin  behavior  as  a 
result  of  variable  thrust,  asymmetric  thrust,  and  thrust 
vectoring . 

Using  a  previously  developed  model  of  the  F-15,  bifurca¬ 
tion  analysis  and  continuation  methods  were  used  to  map  out 
the  equilibrium  and  periodic  solutions  of  the  model  as  a 

function  of  the  thrust  parameters.  A  baseline  bifurcation 
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diagram,  as  a  function  of  alpha  andfjfi^)  of  the  equilibrium 
solutions  for  the  F-15  was  developed.  Thrust  was  varied  and 
changes  were  identified.  Thrust  asymmetries  were  introduced 
and  their  effect  on  entering  and  recovering  from  spins  was 
identified.  Thrust  vectoring  was  introduced  to  see  how 
pitch  and  yaw  vectoring  can  aid  in  the  entry  and  recovery 
from  spins.  Where  deemed  necessary,  time  history  simu¬ 
lations  were  presented  to  further  explain  F-15  behavior. 
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INVESTIGATION  OF  THE  HIGH 
ANGLE  OF  ATTACK  DYNAMICS  OF  THE 
F-15B  USING  BIFURCATION  ANALYSIS 

I .  Introduction 

In  the  military  flight  environment,  pilots  must  regu¬ 
larly  place  themselves  in  the  high  angle  of  attack  flight 
regime  to  out  maneuver  their  opponent.  Unfortunately, 
maintaining  controlled  flight  in  this  regime  is  difficult. 
Loss  of  control  can  occur  through  nonlinear  behavior  such  as 
stalls,  departures,  wing  rock,  nose  slice,  spin  entry,  and 
full  spins.  In  combat  these  will  often  lead  to  fatal 
results.  In  peacetime  training  these  are  not  as  dangerous 
as  long  as  there  is  enough  altitude  to  recover  to  controlled 
flight.  However,  spins  require  quite  a  bit  more  altitude  to 
recover  than  the  other  aircraft  motions.  As  a  consequence, 
many  aircrews  and  their  multimillion  dollar  aircraft  are 
lost  in  spin  accidents.  Between  the  years  1966  and  1970, 
two  hundred  fighter  aircraft  worth  360  million  dollars  were 
lost  in  spin  accidents  resulting  in  100  fatalities  (1:1). 
Even  the  F-15,  the  most  advanced  fighter  in  the  U.S.  Air 
Force  today,  is  not  immune  to  spin  losses.  Baumann  (6:2-4) 
describes  the  details  of  a  recent  Air  Force  accident  inves- 
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tigation  board  in  which  an  F-15  was  lost  due  to  an  flat 
spin.  A  routine  training  flight  turned  into  the  total  loss 
of  an  aircraft  because  of  an  inadvertent  spin. 

SBins 

The  motion  of  an  airplane  in  a  spin  is  characterized  by 
an  angle  of  attack  between  the  stall  and  90  degrees,  and  a 
rapid,  wings  level  descent  toward  the  earth  while  rotating 
about  a  vertical  or  near  vertical  axis  (24:1).  Spins  are 
most  commonly  entered  by  stalling  the  wings  and  introducing 
a  yaw.  The  yaw  increases  the  lift  on  the  wing  outside  of 
the  yaw  and  further  stalls  the  inside  wing.  The  increased 
drag  on  the  stalled  wing  further  drives  the  yaw.  Addition¬ 
ally,  a  rolling  moment  in  the  direction  of  the  yaw  is  intro¬ 
duced  by  the  asymmetric  lift  distribution.  Depending  on  the 
severity  of  the  induced  motions  and  the  aircraft's  physical 
characteristics,  a  spin  may  or  may  not  develop.  A  developed 
spin  is  therefore  a  complex  balance  of  aerodynamic  and  iner¬ 
tia  forces  and  moments.  Once  a  spin  develops,  a  recovery  to 
normal  flight  must  be  accomplished  by  stopping  the  yaw 
rotation  or  breaking  the  stall.  Application  of  a  yawing 
moment  about  the  body  z  axis  opposite  the  spin  is  the  pre¬ 
ferred  method  of  recovery. 


2 


Previous  Studies 


A  rough  idea  of  a  certain  airplane’s  spin  characteris¬ 
tics  can  be  estimated  during  design  by  looking  at  key  aero¬ 
dynamic  and  inertial  factors.  However,  there  is  no  clear 
cut  design  methodology  for  high  angle  of  attack  aerodynamics 
because  of  the  nonlinear  behavior  of  the  fluid  dynamics  of 
separated  flows,  a  high  dependence  on  configuration,  and  a 
lack  of  ground  test  facilities  (10:1).  Vertical  wind  tun¬ 
nels  can  be  used  to  gain  an  idea  of  a  design's  spin  charac¬ 
teristics,  however,  correlation  with  full  scale  aircraft  can 
be  questionable  because  of  Reynolds  number  effects  (24:7). 
Analytical  studies  can  provide  an  additional  tool  in  pre¬ 
dicting  the  spin  characteristics  of  an  airplane. 

Analytical  predictions  of  spin  behavior  have  been 
performed  for  many  years,  in  1954,  Scher  (29)  and  Burk  (8) 
both  demonstrated  that  spins  could  be  simulated  on  computers 
by  using  wind  tunnel  aerodynamic  coefficients  and  solving 
the  nonlinear  equations  of  motion.  Scher  produced  time 
histories  of  spin  entry,  developed  spins,  and  spin  recovery. 
Burk  produced  time  histories  of  spin  recoveries  using  anti¬ 
spin  yaw  moments  and  found  that  the  applied  moment  aided  the 
recovery  from  the  spin.  Additional  time  history  studies 
were  done  in  1959  by  Scher,  Anglin,  and  Lawerence  using  a  60 
degree  delta  wing  airplane  (30),  in  1960  by  Neihouse,  Kli- 
nar,  and  Scher  of  the  X-15  (24),  and  in  1972  by  Adams  of 
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several  airplanes  (1).  In  1966,  Grafton  produced  time 
histories  of  spins  to  determine  the  effect  thrust  has  on 
spins  (13)  and  found  that  generally,  applying  thrust  aided 
the  recovery  from  spins.  These  studies  were  able  to  simu¬ 
late  spins,  but  were  unable  to  accurately  describe  the 
causes  of  nonlinear  behaviors  such  as  jumps  or  the  onset  of 
oscillatory  motion. 

In  1979,  Carrol  and  Mehra  (9)  used  a  uxfferent  approach 
to  analytical  methods  when  they  applied  bifurcation  theory 
and  continuation  methods  to  solve  the  nonlinear  equations  of 
motion.  Equilibrium  solutions  of  the  nonlinear  system  were 
traced  out  by  varying  control  surface  deflections.  The  sta¬ 
bility  was  determined  by  looking  at  the  eigenvalues  of  the 
linearized  system.  The  new  method  was  not  a  simulation, 
like  the  previous  studies,  but  a  map  of  the  airplane's  equi¬ 
librium  solutions  and  therefore  offered  a  more  global  view 
of  the  nonlinear  behavior  of  the  aircraft.  More  important 
though,  the  nature  of  the  transitions  from  stable  to  unsta¬ 
ble  equilibrium  solutions  revealed  the  causes  and  onset  of 
nonlinear  behavior.  Guicheteau  (15,16),  Hui  and  Tobak  (18), 
Zagaynov  and  Goman  (36),  Hawkins  (17),  and  Jahnke  (19,20,21) 
also  used  bifurcation  analysis  to  observe  the  nonlinear 
behavior  of  different  aircraft  configurations,  including  the 
development  of  spins.  Barth  (5)  and  Plansaux  and  Barth  (25) 
investigated  the  nonlinear  behavior  of  the  F-15  using  bifur- 
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cation  analysis  and  proved  its  ability  to  predict  the  air¬ 
craft's  motion,  including  the  onset  of  wing  rock.  However, 
their  model  was  not  realistic  above  40  degrees  angle  of 
attack  and  therefore  could  not  predict  spins.  Beck  (7)  and 
Planeaux,  Beck,  and  Baumann  (26)  continued  this  research 
using  control  augmentation. 

Previous  work  in  P-15  spin  research  using  bifurcation 
theory  was  done  by  Baumann  (6).  He  created  an  F-15  model 
that  more  realistically  describes  the  aero  coefficients  at 
high  angles  of  attack  by  curve  fitting  F-15  aerodynamic  data 
to  angles  of  attack  up  to  90  degrees.  Using  the  new  F-15 
model,  Baumann  found  stable  flat  spins  between  70  and  80 
degrees  angle  of  attack.  These  flat  spins  correlate  well 
with  flight  test  data.  (22,34) 

Overview 

This  paper  continues  the  F-15  spin  research  accom¬ 
plished  by  Baumann  using  a  modified  Baumann  model  which 
includes  control  of  variable  thrust,  asymmetric  thrust,  and 
thrust  vectoring.  Bifurcation  analysis  will  be  used  to 
determine  how  effective  the  thrust  parameters  are  in  causing 
and  recovering  from  flat  spins.  Although  current  opera¬ 
tional  F-15 ' s  do  not  have  the  capability  to  vector  thrust, 
nozzles  with  vectoring  features  are  presently  installed  on 
the  F-15  STOL  demonstrator  aircraft  presently  undergoing 
flight  testing  at  Edwards  AFB,  Ca.  (28:51).  Additionally, 


5 


the  Rockwell /MSB  X-31  aircraft  will  have  thrust  vectoring 
paddles  installed  (27:117).  In  the  X-31  studies,  thrust 
vectoring  will  be  used  to  evaluate  its  combat  utility  at  low 
airspeeds  and  high  angles  of  attack.  Analytical  and  simu¬ 
lation  studies,  such  as  those  done  by  Schneider  (31)  and 
Anderson  (2),  have  shown  the  ability  of  thrust  vectoring  to 
increase  an  aircraft’s  agility.  However,  only  the  study  by 
Burk  (8)  has  gone  on  to  show  the  potential  use  of  applied 
moments  in  spin  recoveries. 

Chapter  11  will  discuss  in  more  details,  the  dynamics  ,o£ 
spins  and  recovery  from  them.  Chapter  III  will  briefly  dis¬ 
cuss  bifurcation  theory  and  the  continuation  method  used  to 
trace  out  the  branches  of  the  bifurcation  diagram.  Chapter 
IV  will  describe  the  F-15  model  and  the  modifications  made 
to  the  Baumann  F-15  model.  Chapter  V  presents  the  results 
that  were  found  during  the  research.  In  Chapter  VI,  the 
conclusions  will  be  presented  and  ideas  of  future  research 
using  this  technique  will  be  given. 
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II. 


Spin  Theory 


The  previous  chapter  defined  spin  motion  and  gave  a  sim¬ 
ple  example  of  how  spins  are  typically  entered.  This  Chap¬ 
ter  will  look  at  flat  spins  and  s-now  briefly  how  the 
aerodynamics  and  physical  characteristics  of  an  airplane 
affect  its  ability  to  recover  from  the  spin.  Most  of  the 
information  in  this  chapter  is  referenced  from  Neihouse, 
Klinar,  and  Scher  (24). 

The  most  dangerous  of  all  spins  is  the  flat  spin.  It  is 
characterised  by  an  angle  of  attack  approaching  90  degrees 
and  high  yaw  rotation  rates.  As  the  plane  approaches  90 
degrees,  the  aerodynamic  control  surfaces  become  ineffective 
due  to  blockage  by  the  wings  and  fuselage  and  recovery  may 
be  difficult.  The  inertia  characteristics  of  modern  fighter 
aircraft  compound  the  difficulty  in  recovering  from  a  flat 
spin.  Since  most  of  the  weight  is  concentrated  in  the  fuse¬ 
lage,  the  high  yaw  rotation  rate  is  accompanied  by  a  large 
amount  of  angular  momentum.  The  aerodynamic  surfaces  must 
provide  moments  to  counter  this  angular  momentum  and  break 
the  spin.  The  decreased  effectiveness  of  the  control  sur¬ 
faces  couples  with  the  large  amount  of  angular  momentum  to 
make  recovery  from  flat  spins  much  more  difficult.  Spin 
test  aircraft  are  often  fitted  with  drogue  parachutes  to  aid 
in  the  recovery  from  flat  spins.  However,  operational  air¬ 
craft  do  not  have  the  luxury  of  spin  parachutes  and  are 


7 


often  lost  due  to  flat  spins.  Vectoring  an  aircraft’s 
thrust  against  the  spin  is  a  potential  source  of  additional 
moments  to  decrease  the  angular  momentum.  Unfortunately, 
modern  fighters  are  not  equipped  with  vectoring  nozzles 
either.  Asymmetric  thrust  settings  in  a  multiengine  air¬ 
craft  are  another  possible  source  for  antispin  moments. 

Since  the  main  component  of  motion  in  a  developed  spin 
is  the  rotation  about  the  z  axis  (yawing  motion),  the  recov¬ 
ery  from  the  spin  therefore  necessitates  decreasing  the  rate 
of  yaw  rotation,  r.  Eqn  (1)  represents  the  r  equation  of 
an  aircraft  in  principal  axes. 


r 


I. 


p  q 


q  S  c 
h 


C» 


p  q 


(l) 


The  yawing  rate  can  be  reduced  by  making  r  negative.  Most 
conventional  aircraft  have  inertia  characteristics  that  make 
the  coupled  term  of  Ix  -  Iy  small.  Therefore,  the  aerody¬ 
namic  yaw  moment  created  by  the  rudder  is  the  primary  moment 
in  recovering  from  spins  for  this  type  of  configuration. 

The  inertia  characteristics  of  modern  fighter  aircraft  make 
the  coupled  term  of  Ix  -  Iy  a  large  negative  number.  A 
decrease  in  the  yaw  rate  can  be  accomplished  by  coupling  the 
roll  and  pitch  rates  to  provide  an  antispin  yaw  moment. 
Assuming  that  the  plane  is  in  a  right  spin,  a  negative  r  ca,' 
be  achieved  by  producing  a  positive  pitch  and  roll  rate.  To 
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most  pilots  this  is  not  intuitive  since  it  requires  prospin 
control  inputs  (i.e.  aileron  into  the  spin  and  pitch  up). 
Accompanying  this  with  an  antispin  rudder  deflection  seems 
to  be  the  most  practical  method  at  recovering  from  a  flat 
spin.  The  F-15  flight  manual  recommends  nearly  full 
aileron/differential  stabilator  deflections  in  the  spin 
direction  to  recover  from  flat  spins  (34:6-7).  The  use  of 
yaw  thrust  vectoring  and  asymmetric  thrust  settings  can  also 
provide  direct  antispin  yaw  moments.  Pitch  vectoring  can 
provide  a  coupling  term  or  can  provide  a  nose  down  moment  to 
try  to  break  the  stall. 
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III.  Bifurcation  Theory 


Nonlinear  phenomena  are  responsible  for  a  variety 
of  effects.  Jumping  between  modes,  sudden  onset  or 
vanishing  of  periodic  oscillations,  loss  or  gain  of  sta¬ 
bility,  buckling  of  frames  and  shells,  ignition,  combus¬ 
tion,  and  chaos  are  but  a  few  examples.  Nonlinear 
phenomena  arise  in  all  fields  of  physics,  chemistry, 
biology,  and  engineering.  The  classical  mathematical 
discipline  that  treats  nonlinear  phenomena  is  bifurca¬ 
tion  theory.  (33, xi) 


This  chapter  will  discuss  the  basic  principles  of  bifur¬ 
cation  theory  and  their  use  in  understanding  nonlinear 
behavior.  Some  of  the  concepts  covered  will  be  equilibrium 
solutions,  stability,  turning  points,  bifurcation  points  and 
Hopf  bifurcations.  Host  of  the  information  in  this  chapter 
is  referenced  from  Seydel  (33).  The  software  program  used 
in  this  research,  AUTO,  will  also  be  described. 


Equilibrium  Points 

Equilibrium  points  represent  steady  states  of  a  dynami¬ 
cal  system;  the  system  is  at  rest  or  in  uniform  motion. 
Equilibrium  points  are  also  referred  to  as  stationary 
points.  The  motion  of  a  non- time  dependent  system  can  be 
modelled  mathematically  as 


u  -  f(a)  (2) 

where  u  is  the  state  vector.  The  equilibrium  states  of  this 
system  would  satisfy  the  equation 
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f (u)  -  o. 


(3) 


An  aircraft  in  this  state  would  exhibit  no  translational  or 
angular  accelerations  and  would  have  constant  roll  and  pitch 
angles . 

The  dependence  of  the  system  on  some  control  parameter 
can  be  found  by  varying  the  parameter  and  finding  any  new 
equilibrium  points.  The  equation  can  be  represented  as 

-f(u.X)  (4) 

where  \  is  the  control  parameter.  This  parameter  is  called 
the  bifurcation  parameter.  In  an  aircraft  model,  these 
parameters  would  be  such  things  as  the  elevator  deflection 
or  thrust  level.  A  qualitative  idea  of  system's  dependence 
on  the  varying  parameter  is  found  by  plotting  the  new  value 
of  a  representative  state  variable  versus  the  value  of  the 
parameter.  This  diagram  is  called  a  bifurcation  diagram. 

For  the  aircraft  model  with  the  elevator  deflection  as  the 
bifurcation  parameter,  a  bifurcation  diagram  provides  an 
idea  of  the  aircraft's  equilibrium  motion  as  the  elevator  is 
deflected  from  one  value  to  another.  If  the  elevator  is 
varied  from  stop  to  stop,  the  global  behavior  of  the 
aircraft  can  be  found.  Unfortunately,  these  diagrams  tell 
little  of  the  aircraft  response  over  time  or  of  the. 
stability  of  the  equilibrium  points. 
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Stability 


The  stability  of  an  equilibrium  point  is  determined  by 
identifying  whether  the  system  will  return  to  the  equilib¬ 
rium  point  if  it  is  disturbed.  A  point  is  considered  stable 
if  the  response  to  a  small  perturbation  is  small  as  time 
goes  to  infinity.  If  the  response  goes  to  zero  as  time  goes 
to  infinity  it  is  considered  asymptotically  stable.  It  is 
unstable  if  the  response  grows  as  time  goes  to  infinity.  A 
neutrally  stable  point  would  neither  go  to  zero  nor  grow. 

The  size  of  the  perturbation  is  important  since  an  equilib¬ 
rium  point  may  be  stable  for  a  small  perturbation  but  unsta¬ 
ble  for  a  larger  one. 

Stability  can  be  found  by  linearizing  the  system  around 
the  equilibrium  point.  This  is  a  good  approximation  of  the 
nonlinear  system  close  to  the  equilibrium  point.  The  sta¬ 
bility  can  then  be  determined  by  looking  at  the  eigenvalues 
of  the  Jacobian  matrix  of  the  linear  system.  The  system  is 
stable  if  the  real  parts  of  the  eigenvalues  are  negative  or 
zero.  A  positive  real  part  indicates  that  the  point  is 
unstable. 

Turning  Points 

A  turning  point  in  a  nonlinear  system  has  a  single 
eigenvalue  equal  to  zero.  Figure  3-1  is  a  bifurcation  dia¬ 
gram  with  a  turning  point  found  in  the  differential  equation 
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y 


(5) 


~  \  —  y 2 . 

At  the  turning  point  or  limit  point,  the  only  equilibrium 
solution  is  X.  =  0,  y  =  0,  With  \  >  0  there  are  two 
solutions.  The  solution  +  is  stable,  while  -fh  is 
unstable.  Turning  points  do  not  necessarily  separate  stable 
equilibria  from  unstable  equilibria.  Unstable  solutions  can 
also  exist  on  both  sides  of  the  turning  point  as  one 
eigenvalue  crosses  zero.  Fig.  3-1  also  identifies  stable 
equilibria  as  being  solid  lines  while  unstable  equilibria 
are  identified  by  dashed  lines. 


Figure  3-1  Turning  Point  on  a  Bifurcation  Diagram 

A  unique  phenomenon  called  hysteresis  occurs  when  a 
branch  loses  stability  at  a  turning  point  and  then  becomes 
stable  at  another  turning  point.  Fig.  3-2  is  a  bifurcation 
diagram  of  a  typical  hysteresis  point.  Hysteresis  leads  to 
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jump  phenomena  between  the  stable  branches  as  the  parameter 
is  varied  beyond  either  limit  point.  Aircraft  can  display 
jump  phenomena  between  stable  states,  such  as  the  jump  from 
low  angle  of  attack  equilibria  to  a  high  angle  of  attack 
spin.  Two  good  articles  dealing  with  jump  phenomena  in  air¬ 
craft  maneuvers  are  written  by  Schy  and  Hannah  (32)  and 
Young,  Schy,  and  Johnson  (35). 


X 

Figure  3-2  Limit  Points  Showing  Hysteresis 


Bifurcation  Points 

Bifurcation  points  also  have  one  zero  eigenvalue.  How¬ 
ever,  unlike  a  turning  point,  there  are  solutions  for  values 
of  X  on  both  sides  of  a  bifurcation  point.  Fig.  3-3  is  an 
example  of  the  pitchfork  bifurcation  that  develops  in  the 
differential  equation 

y  -  Xy  -  y3.  (6) 
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For  all  values  of  X.  the  trivial  solution,  y  =  0,  is  an 
equilibrium  solution.  Additionally,  for  \  >  0  there  are  two 
nontrival  equilibria,  y  =  ±Jx.  However,  the  branch  y  =  0 
losses  stability  at  the  bifurcation  point  and  a  bifurcation 
of  two  stable  branches  occurs.  This  is  referred  to  as  a 
supercritical  pitchfork. 


If  the  branch  y  =  0  gains  stability  at  the  bifurcation  point 
and  a  bifurcation  of  unstable  branches  occurs,  the  result  is 
called  a  suhcritical  pitchfork.  Fig.  3-4  is  an  example  of  a 
subcritical  bifurcation.  The  behavior  at  these  bifurcation 
points  is  also  referred  to  as  a  stationary  or  steady  state 
bifurcation. 
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Figure  3-4  Subcritical  Pitchfork 


Hopf  Bifurcation 

The  types  of  nonlinear  phenomena  identified  so  far  are 
for  equilibrium  solutions.  The  type  of  bifurcation  that 
connects  equilibrium  solutions  with  periodic  motion  is  the 
Hopf  bifurcation.  Periodic  solutions  arise  at  points  where 
two  eigenvalues  of  the  linearized  system  become  purely  imag¬ 
inary.  For  an  example,  a  Hopf  bifurcation  arises  from  the 
two  equations 


y i  -  -Yz  +  yiC^  -  y \  -  yi) 

Yz  -  Yi  +  y2(^  -  y\  -  yi)- 

The  only  equilibrium  solution  for  all  K  occurs  at 
y i  -  y2  -  0.  However,  the  eigenvalues  are  which 
indicates  that  the  equilibrium  points  are  unstable  for  K  >  0 
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and  stable  for  X  <  0.  A  Hopf  point  is  located  at  X  =  0 
since  both  eigenvalues  are  purely  imaginary.  Additionally, 
a  r’-ange  of  stability  takes  place  without  a  turning  point  or 
any  branches  bifurcating.  The  exchange  of  stability  occurs 
through  the  formation  of  a  family  of  limit  cycles  at  the 
Hopf  point. 

Limit  cycles  are  found  by  putting  yi  and  y2  into  polar 
coordinates 


y,  -  pcosO,  y2  -  psin0  (8) 

and  then  substituting  them  into  Eqn  (7).  By  manipulation 
this  yields 


p  -  p(\  -  p2)  (9) 

e  -  i.  (io) 

This  shows  that  6=1  which  is  not  an  equilibrium  solution. 
For  X.  >  0  the  result  is  a  periodic  orbit  with  an  amplitude 
growing  by  %/x.  Fig.  3-5  shows  how  this  looks  in  three 
dimensions.  The  limit  cycle  encircles  the  unstable  equilib¬ 
rium.  Fig.  3-6  is  an  example  of  a  Hopf  bifurcation  on  a 
bifurcation  diagram.  A  stable  branch  goes  unstable  at  the 
Hopf  point.  The  circles  represent  the  maximum  amplitude  of 
the  limit  cycle  at  X.  Closed  circles  indicate  stable  limit 
cycles  and  open  circles  indicate  unstable  limit  cycles. 

Periodic  solutions  lose  stability  via  three  mechanisms; 
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Figure  3-5  Limit 


turning  points,  period  doubling,  and  bifurcation  into  a 
torus.  Floquet  multipliers,  which  are  analogous  to  eigenva¬ 
lues,  are  used  to  find  the  stability  of  a  limit  cycle. 

AUTO  Software 

The  tool  used  in  this  research  to  trace  out  equilibrium 
branches,  determine  stability,  identify  turning  points  and 
bifurcation  points,  and  find  limit  cycles  is  the  program 
AUTO  written  by  Doedel  (11).  From  a  known  starting  point, 
(u0,X.0),  which  satisfies  Eqn  (3),  Doedel  uses  the  psuedo 
arclength  continuation  technique  to  trace  out  equilibrium 
solutions  for  new  values  of  X..  The  psuedo  arclength  tech¬ 
nique  varies  the  stepsize  along  the  branch  and  using  the 
direction  vector  (u.k)  a  predictor-corrector  algorithm 
finds  the  next  solution.  The  predictor/corrector  algorithm 
used  is  the  Newton  method.  The  psuedo  arclength  technique 
allows  the  algorithm  to  be  scaled  so  it  can  compute  near  and 
past  limit  points  where  the  direction  vector  is  infinite. 
Doedel  also  incorporates  an  adaptive  stepsize.  If  the  solu¬ 
tion  converges  rapidly  using  the  predictor/corrector  algo¬ 
rithm,  the  stepsize  is  increased  to  save  computation  time. 
Additionally,  if  the  solution  does  not  converge,  the 
stepsize  is  halved  until  a  minimum  stepsize  is  reached.  The 
program  will  then  signal  nonconvergence. 

AUTO  identifies  bifurcation  points  and  turning  points  by 
monitoring  the  Jacobian  matrix  at  each  solution  and  identi- 
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fying  sign  changes  in  the  eigenvalues.  Using  bifurcation 
analysis,  AUTO  identifies  these  changes  as  limit  points, 
bifurcation  points,  or  Kopf  Bifurcations.  AUTO  continues  on 
the  main  branch  until  a  user  specified  number  of  points  is 
reached  or  values  of  K  or  u  exceed  user  specified  limits. 
AUTO  has  the  capability  to  go  back  to  the  bifurcation  points 
to  compute  the  branches  emanating  from  the  bifurcation 
point.  Additionally,  AUTO  can  go  back  and  compute  the  limit 
cycles  that  begin  at  the  Hopf  bifurcations.  More  informa¬ 
tion  on  the  capabilities  of  AUTO  can  be  found  in  the  AUTO 
user  manual  (11). 
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IV.  Model  Development 


Aircraft  Description 

The  P-15B  aircraft,  which  is  modelled  in  this  research, 
is  a  two  seat,  high  performance,  supersonic,  all-weather 
air-superiority  fighter.  The  aircraft's  primary  mission  is 
aerial  combat,  however,  it  can  also  be  configured  for  ground 
attack.  It  is  powered  by  twin  Pratt  and  Whitney  F-100  tur¬ 
bofan  engines.  The  model  developed  for  this  research 
includes  pitch  and  yaw  thrust  vectoring  which  a  baseline 
F-15B  does  not  have.  The  vectoring  nozzles  are  assumed  to 
have  no  effects  on  the  aerodynamic  characteristics  or  weight 
and  balance  of  the  F-15B.  Appendix  A  provides  the  physical 
dimensions  of  the  aircraft  and  weight  and  balance. 

The  aircraft's  aerodynamic  control  surfaces  are  the 
ailerons,  rudder,  elevator.  Thrust  settings  can  be  indepen¬ 
dently  controlled  for  both  right  and  left  engines.  For  this 
research,  the  yaw  and  pitch  angles  of  the  nozzles  can  also 
be  controlled.  Several  control  characteristics  of  the  F-15B 
will  not  be  modelled  to  simplify  the  research.  These 
include  the  effects  of  the  Control  Augmentation  System 
(CAS),  Aileron  Rudder  Interconnect  (ARI ) ,  and  speedbrake. 

The  differential  elevator  deflection  is  also  set  at  a  con¬ 
stant  gain  times  the  aileron  deflection.  The  F-15  aero 
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coefficients  modelled  are  for  low  speed  flight  and  constant 
altitude.  Therefore,  flight  conditions  of  20,000  feet  and 
low  Mach  numbers  will  be  used. 

Force  and  Moment  Equations 

The  force  and  moment  equations  used  in  this  research 
are  the  body-axis  force  and  moment  equations  used  by  Baumann 
(6:20-21)  but  modified  to  include  forces  and  moments  due  to 
variable  thrust,  asymmetric  thrust  and  thrust  vectoring. 
These  equations  therefore  have  both  aerodynamic  and  thrust 
components.  The  aerodynamic  force  and  moment  coefficients 
were  modelled  for  the  F-15  by  Baumann.  He  used  a  statis¬ 
tical  software  program  to  curve  fit  F-15  aerodynamic  data 
from  -20  degrees  to  90  degrees  angle  of  attack.  These  curve 
fits  are  fairly  representative  of  the  actual  F-15  aerody¬ 
namic  coefficients  (23).  Some  amount  of  data  smoothing  can 
be  expected,  however,  the  general  trends  in  the  data  are 
maintained.  The  details  of  his  curve  fitting  techniques  can 
be  found  in  (6).  At  low  angles  of  attack,  the  above  coeffi¬ 
cients  are  symmetric  with  respect  to  the  lateral  variables 
(3,  p,  and  r.  However,  asymmetries  are  present  above  40 
degrees  angle  of  attack  due  to  asymmetric  shedding  of  nose 
vortices  (22:3.4). 

The  thrust  components  are  added  to  the  aerodynamic  coef¬ 
ficients  to  produce  combined  thrust/aerodynamic  force  and 
moment  coefficients.  The  thrust  contributions  to  the 
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modified  force  coefficients  are  found  by  determining  the 
thrust  in  the  body  x,y,  and  z  directions  as  a  result  of 
total  thrust  and  any  vectoring.  The  moments  are  found  by 
determining  the  contribution  by  each  engine  in  the  body  x,  y 
and  z  direction  and  multiplying  each  by  the  appropriate 
moment  arm  to  the  center  of  gravity.  Asymmetric  thrust  is 
modelled  using  the  variables  right  engine  thrust  (Tr)  and 
left  engine  thrust  (T,) .  Pitch  vectoring  is  modelled  using 
the  variable  6pv  with  a  positive  value  causing  the  a  nose  up 
pitching  moment  (flow  deflected  upward).  Yaw  vectoring  is 
modelled  using  the  variable  6yv  with  a  positive  value  caus¬ 
ing  a  yaw  to  the  left  (flow  deflected  left).  The  moment  arm 
offsets  are  defined  as  dTx.  dTy.  and  dT».  Fig.  4-1  shows 
these  variables. 


Figure  4-1  Physical  Description  of  Thrust  Variables 
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Etkin  (12)  was  used  to  verify  Baumann's  equations  and  the 
thrust  effects.  The  resulting  equations  are: 


Cx  “  CL(a,6,)  sina  -  C0(a,6,)  cosa  +  Tx/(q  S)  (11) 

Cy  =  Cy(0t||B|,6,)  +  Cy8,(a)  6,  +  Cy5r(<l)6r 

+  [b/2Vtr)[CyrCa)r  +  C„(a)p]  +  ACy>.(a,  P)  (12) 

+  Cy#^,(ot ,  6 , )  6 +  Ty/(q,  S) 

Cx  -  -  CL(a,6,)  cosa  -  C0(a,6,)  sina  +  Tt/(q  S)  (13) 

Ci  “  C1(J(a,  |  p  |)  p  +  C18,(a,  6#)  6,  +  C)5r(a,|6r|)  6r 

+  [b/2Vtr][Clp(a)  p  +  Clr(a)  r]  (14) 

+  Clta,(a,6.)  6a,  +  AC|^(a,P) 

+  [T„dTy  -  TEldTy]/(q  S  b) 

cm  -  CB0(a,6.)  +  [c/(2  Vtt)JC-q(a)  q  (15) 

+  [T xdTl -  T *dTx]/ (q  S  c) 


C„  =  Cn,(a.  IPI.6.)  p  +  Cn#,(a)  6,  +  Cu0r(a.  |  p  1. 1  6,  1.6.)  6r 

+  [b/2V&][Cnp(a)  p  +  CBr(a)  r]  +  CBBBe(a,  6Ae)6Ae  (16) 

+  ACD,(a.P)  +  ACni.(a,  P) 

+  [1yd  T*  ~  1  xr d  Tr  “  Txldrr)^1!  S  b). 

These  coefficients  are  used  in  the  equations  of  mo. ion. 
Equations  of  Motion 

The  equations  of  motion  for  an  airplane  are  r’eiived 
from  Newton's  Law 


Ma 

(17) 

la 

(18) 
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where  a  is  the  aircraft  translational  acceleration,  a  is  the 

aircraft  rotational  acceleration,  M  is  the  aircraft  mass,  I 
is  the  aircraft  rotational  inertia  tensor,  and  F  and  N  are 
the  forces  and  moments  acting  on  the  aircraft.  The  result 
is  a  twelfth  order  system.  However,  by  making  the  following 
assumptions,  rigid  aircraft,  constant  air  density,  constant 
gravity,  and  a  flat  earth,  and  by  transforming  the  force  and 
moment  equations  from  the  inertial  frame  to  a  frame  fixed  to 
the  aircraft,  the  equations  of  motion  are  reduced  from  order 
twelve  to  order  nine.  The  aircraft  state  can  be  described 
by  the  nine  state  variables  (a,  p,  p.  q,  r,  0.  4>.  tj>,  V).  If 
the  xz  plane  is  a  plane  of  symmetry,  the  following  equations 
are  formed.  Translational  acceleration  equations: 
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Rotational  acceleration  equations: 
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The  aircraft  orientation  or  Euler  angle  equations  are: 

0  -  qcos#  -  rsin$ 

4»  -  P  *■  (qsin$  +  rcos$)tan0 
V  -  Cqsin4>  +  rcos#)sec0. 
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The  yaw  angle  is  decoupled  from  the  rest  of  the  equations 
via  the  definition  of  the  Euler  angles.  The  aircraft  is 
first  rotated  through  the  yaw  axis,  then  the  pitch  axis  and 
then  the  roll  axis.  The  yaw  rotation  does  not  change  the 
direction  of  the  gravity  vector  since  the  z  axis  of  the 
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aircraft  and  earth  are  initially  aligned.  The  result  is  an 
eighth  order  model  without  the  orientation  equation.  The 
above  equations  are  in  the  form  of  Eqn  (4)  with  the  state 
vector  u  -  [a,  p,  p,  q,  r,  0,  V]T  and  the  variable  param¬ 
eters  \  -  [5,,  6r,  8,.  Tl,  Tr,  6pv,  8rv]T.  The  purpose  of  this 
research  is  to  find  the  solutions  to  this  set  of  equations 
that  satisfy  Eqn  (3). 

Model  Modifications 

In  addition  to  the  modifications  due  to  the  thrust 
variables,  one  curve  fit  in  the  Baumann  model  was  revised. 
Initial  elevator  sweeps  of  the  Baumann  model  identified  a 
loss  of  longitudinal  stability  at  low  angles  of  attack  and 
then  a  regain  of  stability  shortly  later.  This  loss  of 
stability  did  not  match  flight  test  data.  The  curve  fits 
for  the  longitudinal  aerodynamic  coefficients  were  compared 
to  the  aero  database  to  see  if  they  were  in  error.  The 
pitch  damping  derivative  was  identified  as  the  cause  of  the 
problem.  The  curve  fit  of  (!„„  in  Baumann’s  model  was  inac¬ 
curate  at  both  low  angles  of  attack  and  above  70  degrees.  A 
new  curve  fit  was  developed  to  replace  the  incorrect 
equation  for  Coq.  Pig.  4-2  compares  the  original  and  modi¬ 
fied  curve  fits  for  Cmq  and  also  provides  an  idea  as  to  how 
the  aero  data  is  modelled. 
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V.  Results 


The  main  point  of  this  research  was  to  see  the  effects 
of  thrust,  thrust  vectoring,  and  asymmetric  thrust  on  the 
spin  characteristics  and  spin  recovery  of  the  F-15.  There¬ 
fore,  a  baseline  of  the  global  spin  characteristics  of  the 
F-15  had  to  be  established  to  compare  the  results  against. 
The  Baumann  model  was  initially  run,  however  some  irregula¬ 
rities  in  the  data  at  low  alpha  necessitated  changes  in  the 
model .  The  model  used  in  this  research  is  a  modified 
Baumann  model .  The  new  model  is  first  compared  to  the  old 
model  to  see  the  effects  of  the  changes.  The  new  model  will 
then  be  further  evaluated  to  define  the  global  F-15  spin 
characteristics.  Thrust  levels  will  be  varied  to  see  how 
they  affect  these  spin  characteristics.  Asymmetric  thrust 
will  be  introduced  to  see  if  it  can  lead  to  spins  and  to  see 
if  it  can  be  used  to  recover  from  flat  spins.  Finally, 
pitch  and  yaw  thrust  vectoring  will  be  used  to  see  how  they 
can  aid  in  the  recovery  from  flat  spins. 

The  general  technique  used  in  the  investigation  was  to 
find  the  global  effects  of  varying  the  parameter  of  inter¬ 
est,  and  then  to  concentrate  on  the  stable  flat  spin  regime. 
Additional  bifurcation  diagrams  were  developed  from  a 
starting  point  in  the  stable  flat  spin  region.  Time  histo¬ 
ries  are  shown  in  some  cases  to  provide  a  better  understand¬ 
ing  of  the  F-15's  behavior. 
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Comparison  with  Unmodified  Model 


The  first  run  in  this  research  was  an  elevator  sweep 
from  a  starting  point  calculated  by  Baumann  (6:32-34).  Bau¬ 
mann  used  a  thrust  setting  of  8300  lbf  in  his  research 
because  it  is  the  thrust  setting  for  trim  conditions  at  0.6 
Mach  and  20,000  feet  (6:19).  The  a-6,  bifurcation  diagram 
from  this  elevator  sweep  indicated  a  small  area  of  instabil¬ 
ity  bounded  by  Hopf  points  at  a  fairly  low  angle  of  attack. 
As  discussed  earlier,  an  incorrect  curve  fit  for  Cmq  was 
identified  as  the  cause  of  the  instability  and  also  showed 
that  the  original  curve  fit  did  not  represent  the  high  angle 
of  attack  data.  A  new  curve  fit  was  found  for  Cmq  and  was 
incorporated  into  the  model.  Again  an  elevator  sweep  was 
run  from  the  original  starting  point.  The  a-6,  bifurcation 
diagram  for  the  modified  model  did  not  have  the  small  unsta¬ 
ble  area.  A  comparison  of  the  Baumann  and  modified  model  in 
the  high  a  regime  was  then  accomplished  using  a  rudder  sweep 
from  a  starting  point  near  full  elevator  deflection.  The 
original  elevator  sweep  did  not  continue  into  the  spin 
region  and  therefore  Baumann  found  that  rudder  sweeps  can  be 
used  to  reach  this  area  (6:47).  Since  the  new  curve  for  C„,q 
is  not  as  stable  in  pitch  damping  at  high  angles  of  attack 
as  the  old  model  (see  Fig.  4-2),  the  new  model  should  show  a 
smaller  region  of  stability  in  the  spin  region.  The  a-6r 
bifurcation  diagram  of  the  unmodified  model  (Fig.  5-1)  shows 
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a  small  region  of  stable  spins  with  full  elevator  deflec¬ 
tion.  The  a  -  6r  diagram  using  the  new  model  (Fig.  5-2)  shows 
the  same  tracing  of  equilibrium  solutions  as  Fig.  5-1. 
However,  the  stable  portion  of  the  branch  is  no  longer  pres¬ 
ent.  This  matches  the  expected  results  that  the  equilibrium 
would  be  less  stable.  A  stable  limit-cycle  may  exist  in  the 
modified  model  in  the  vicinity  of  the  stable  area  of  the 
unmodified  model.  However,  the  rudder  sweep  was  not  able  to 
find  these  limit  cycles  because  there  are  no  Hopf  points 
present.  With  full  elevator  deflection,  the  new  model  is 
similar  to  the  unmodified  model  except  for  the  loss  of  the 
stable  spin  branch. 
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Figure  5-1  Baumann  Model  Rudder  Sweep 
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Figure  5-2  Revised  Model  Rudder  Sweep 

Basel ine  Model  Spin  Characteristics 

The  previous  comparison  was  done  by  holding  the  elevator 
and  ailerons  constant  and  varying  the  rudder.  This  intro¬ 
duces  asymmetries  into  the  airplane  aerodynamics  by  varying 
the  rudder.  A  different  method  of  producing  bifurcation 
diagrams  is  to  hold  neutral  ruddar  and  aileron  and  vary  the 
elevator.  This  method  provides  a  global  view  of  the  F-15 
longitudinal  motion  in  a  symmetric  configuration  as  a  func¬ 
tion  of  elevator  deflection.  Any  asymmetries  identified 
would  therefore  be  the  result  of  aerodynamic  asymmetries  and 
not  the  result  of  rudder  or  aileron  inputs.  Rudder  sweeps 
are  still  necessary,  however,  to  provide  starting  points  for 
high  angle  of  attack  branches.  Fig.  5-3  is  the  baseline 
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Figure  5-2  Revised  Model  Rudder  Sweep 
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a-6,  bifurcation  diagram  of  the  P-15B  for  a  thrust  level  of 
8300  Ibf.  This  diagram  shows  only  the  equilibrium  solutions 
found  between  the  F-15  elevator  deflection  limits.  The 
entire  continuation  diagram  can  be  found  in  Appendix  C,  Fig. 
C-l. 

Several  equilibrium  branches  are  identified  in  Fig.  5-3. 
The  low  angle  of  attack  stable  branch  loses  stability  at  6. 

=  -19.3  degrees  and  identifies  the  onset  of  wing  rock.  The 
unstable  branches  at  alpha  =  36  degrees  and  between  alpha  = 
40  degrees  and  50  degrees  identify  spirals.  The  branch 
found  between  alpha  =  64  degrees  and  84  degrees  is  the  spin 
branch.  It  contains  both  stable  and  unstable  equilibrium 
and  limit-cycle  solutions.  The  stable  portion  of  the  branch 
between  6.  =  -21.4  degrees  and  -15.4  degrees  identifies  a 
stable  right  spin.  It  is  bounded  by  stable  limit-cycles 
that  eventually  become  unstable.  The  unstable  equilibrium 
branch  connected  to  the  stable  branch  has  a  maximum  value  of 
the  unstable  eigenvalue  equal  to  0.10683  =*>  2.4174i  at  6,  =  0 
degrees.  This  is  important  since  the  positive  eigenvalue  is 
a  small  number.  This  means  that  the  loss  of  stability  along 
the  branch  will  most  likely  not  be  an  immediate  event  (i.e. 
long  transient).  The  small  stable  region  at  6.  =  -25 
degrees  identifies  a  stable  left  spin.  It  gains  stability 
via  a  turning  point  and  loses  stability  via  a  Hopf  bifurca¬ 
tion  and  with  stable  limit-cycles.  The  larger  area  of  right 
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spin  with  symmetric  controls  can  be  explained  by  the  asym¬ 
metric  yawing  moment  above  alpha  =  4C  degrees  (22:3.4).  The 
asymmetry  favors  a  right  yaw  and  therefore,  a  larger  right 
spin  region.  Of  additional  interest  is  that  the  stable  spin 
branch  lies  directly  above  point  were  wing  rock  begins.  The 
F-15  could  theoretically  encounter  wing  rock  and  then  jump 
to  the  stable  spin  branch. 

The  F-15  flight  manual  recommends  several  methods  to 
recover  from  spins  (34:6-7).  The  highly  oscillatory  spin 
can  be  recovered  by  neutralizing  the  controls.  Stable  flat 
spins  can  be  recovered  by  applying  aileron  in  the  direction 
of  the  spin.  As  discussed  earlier,  by  applying  aileron  in 
the  direction  of  the  spin,  a  cross  coupling  inertia  effect 
acts  in  the  direction  opposite  the  yaw.  The  manual  also 
states  that  rudder  deflection  in  either  direction  has  little 
effect  on  spin  recovery.  This  was  also  discussed  earlier 
and  is  a  result  of  the  rudder  being  washed  out  by  the  wake 
off  the  wings  and  fuselage.  These  recovery  techniques  were 
applied  in  a  simulation  of  the  stable  flat  spin  at  6.  = 

-19.14  degrees  to  see  if  the  model  corresponds  to  actual 
flight  behavior.  The  selection  of  6,  =  -19.14  degrees  as  a 
starting  point  is  motivated  by  the  fact  that  it  is  near  the 
center  of  the  stable  branch  and  corresponds  to  the  onset  of 
wing  rock  at  low  a.  This  starting  point  will  be  used  in 
many  of  the  bifurcation  diagrams  in  the  following  sections. 
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Fig.  5-4  shows  the  a  -  time  simulation  using  elevator  for 
recovery.  If  the  elevator  is  neutralized,  the  F-15  will 
enter  an  oscillatory  spin  at  a  lower  angle  of  attack.  If 
the  elevator  is  deflected  to  a  negative  full  deflection,  the 
F-15  enters  an  oscillatory  spin  at  a  higher  angle  of  attack. 
The  F-15  will  not  recover  from  the  flat  spin  in  a  reasonable 
time  using  elevator  alone.  Comparing  Fig.  5-4  to  Fig.  5-3, 
the  oscillatory  spins  begin  at  the  Hopf  points  bounding  the 
stable  region.  These  oscillations  are  centered  on  the 
unstable  equilibrium  branches  and  can  be  classified  as 
unstable  oscillations  because  they  continue  to  grow. 
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Pig.  5-5  is  the  simulation  using  rudder  for  recovery.  A 
fully  deflected  antispin  rudder  recovers  the  F-15  to  wing 
rock  in  13  seconds.  A  fully  deflected  prospin  rudder  drives 
the  F-15  into  an  unstable  oscillatory  spin  that  recovers  the 
F-15  to  wing  rock  in  55  seconds.  In  a  real  world  situation 
the  F-15  would  not  have  55  seconds  to  recover  from  a  flat 
spin,  however  this  diagram  is  valuable  in  revealing  the  true 
nature  of  the  oscillatory  spin  as  being  unstable.  The 
effectiveness  of  the  rudder  to  recover  from  a  flat  spin  in 
this  simulation  is  contrary  to  the  statement  in  the  F-15 
flight  manual . 


Figure  5-5  Simulation  of  Spin  Recovery  Using  Rudder 
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Fig.  5-6  is  the  recovery  attempt  using  ailerons-  Aileron  in 
the  direction  of  the  spin  actua.ly  increase  the  angle  of 
attack  of  the  spin  while  aileron  opposite  the  spin  reduces 
the  angle  of  attack.  The  F-15  should  not  be  able  to  recover 
using  ailerons  alone.  The  results  of  this  simulation  are 
also  contrary  to  the  statements  in  the  F-15  flight  manual, 
at  least  for  6.  =  -19.14  degrees  and  a  thrust  level  of  830C 
lbs . 


Figure  5-6  Simulation  of  Spin  Recovery  Using  Ailerons 

These  discrepancies  were  investigated  by  looking  at  the 
equations  for  the  aero  coefficients  due  to  the  aileron  and 
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rudder  deflections  and  seeing  if  they  were  the  cause.  The 
equations  of  motion  were  double  checked  to  see  if  they  were 
correct.  Possible  cross  coupling  effects  were  investigated 
by  looking  at  the  stable  spin  angular  velocities. 

The  rudder  effectiveness  is  caused  by  the  curve  fit  for 
Cn4r.  The  curve  fit  is  dependent  on  6r,  a,  p,  and  6,  with 
generally  any  value  of  (3  and  negative  values  of  6,  increasing 
the  rudder  effectiveness  at  high  a.  Unfortunately,  the 
curve  fit  makes  the  rudder  deflections  twice  as  effective  at 
large  negative  values  of  6,  than  aerodynamic  data  indicates. 
Therefore,  the  model  is  somewhat  inaccurate  in  the  effects 
of  rudder  at  high  angles  of  attack.  The  starting  points 
used  for  the  continuation  of  the  high  angle  of  attack 
branches  for  the  a  -  6,  bifurcation  diagram  are  still  accu¬ 
rate  since  at  these  points  the  rudder  in  undeflected. 

The  ineffectiveness  of  the  ailerons  can  be  explained  as 
a  cross  coupling  effect  in  the  model.  Table  1  shows  the 
values  of  p,  q,  and  r  in  the  stable  spin  regime  for  6,  = 
-19.14  degrees. 

Table  I.  Stable  Spin  Angular  Velocities 


p  rad/sec 

q  rad/sec 

r  rad/sec 

0.6119 

-0.0874 

1.932 
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Since  I*  -  Iy  for  the  F-15  is  negative,  it  couples  with  a 
positive  p,  and  negative  q  to  drive  the  yaw  rate  -  see  Eqn 
(1).  Therefore,  in  the  right  flat  spin  region  being 
investigated,  prospin  aileron  deflections  are  not  effective 
in  recovering  from  the  spin  and  actually  drive  the  yaw  rate. 
Deflections  opposite  the  spin  can  slow  the  spin  rate, 
however,  there  is  not  enough  control  authority  to  completely 
stop  it.  This  anomaly  is  therefore  a  consequence  of  the 
spin  characteristics  and  not  a  problem  with  the  model. 

Since  the  main  investigation  in  this  research  is  thrust 
effects,  the  discrepancy  due  to  rudder  and  aileron  is  noted 
but  will  not  be  further  investigated.  The  model  will  be 
kept  in  a  symmetric  6r,  6.  configuration  during  investigation 
of  the  thrust  effects  and  therefore  these  discrepancies  will 
not  enter  into  the  rest  of  the  results.  The  coupling 
effects  will  be  important  during  pitch  vectoring,  and  will 
be  discussed  more  in  that  section. 

Throttling 

The  effects  of  varying  the  thrust  level  in  the  model  are 
discussed  next.  First,  throttling  was  run  from  the  origi¬ 
nal  level  of  8300  lbf  to  see  how  the  model  reacted  at  low 
angles  of  attack.  The  elevator  deflection  was  set  at  -19.14 
degrees  with  neutral  rudder  and  ailerons.  The  a  -  T  bifur¬ 
cation  diagram  is  shown  in  Fig.  5-7.  The  diagram  is  inter¬ 
esting  in  that  a  does  not  change  with  thrust  until  the 
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thrust  to  weight  ratio  approaches  one.  At  37,000  lbs 
thrust,  the  diagram  hits  a  wall  and  shoots  to  a  =  90 
degrees.  For  thrust  to  weight  ratios  greater  than  1,  the 
model  is  in  a  constant  acceleration  and  there  could  be  no 
solutions  that  would  satisfy  Eqn  (3).  From  this  diagram, 
however,  two  starting  points  were  picked  to  continue  up  to 
high  angles  of  attack  using  rudder  sweeps. 
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The  two  values  selected  for  continuation  were  0  Ibf  and 
29200  lbf.  The  zero  thrust  value  was  picked  to  show  a  base¬ 
line  of  the  spin  region  with  no  thrust.  The  thrust  value  of 
29200  lbf  was  selected  since  it  is  close  to  the  full 
military  power  rating  (non-afterburning)  of  the  PW  F-100 
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Figure  5-7  Low  Alpha  Engine  Throttling 
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engines.  Rudder  sweeps  were  again  used  to  provide  high 
angle  of  attack  starting  points  for  the  elevator  sweeps. 

The  longitudinal  motion,  symmetric  spin  characteristics  were 
then  found  by  doing  elevator  sweeps  with  the  thrust  held  at 
0  lbf  '.Fig .  5-8)  and  at  29200  lbf  (Fig.  5-9).  The  full 
continuation  diagrams  of  these  elevator  sweeps  can  be  found 
in  Appendix  C. 

Comparing  the  three  elevator  sweeps  for  different  thrust 
levels  (Figs.  5-3,  5-8,  and  5-9)  the  stable  right  spin 
branch  decreases  in  size  with  thrust.  This  result  can  be 
partly  attributed  to  a  nose  up  moment  due  to  thrust.  The 
engines  are  located  below  the  aircraft’s  center  of  gravity 
and  therefore  a  pitch-up  moment  accompanies  increases  in 
thrust.  This  increases  the  value  of  q  and  the  cross  cou¬ 
pling  effect  of  q  and  r  would  bring  about  an  antispin 
moment.  A  better  understanding  of  this  phenomenon  is 
achieved  by  creating  a  bifurcation  diagram  of  a  -  T  from  a 
starting  point  in  the  stable  spin  region.  Fig.  5-10  shows 
that  thrust  decreases  the  angle  of  attack  and  eventually 
produces  limit  cycle  behavior  with  high  thrust  levels,  how¬ 
ever  it  will  not  bring  the  plane  out  of  the  spin.  A  simu¬ 
lation  was  run  to  show  this  effect  and  is  shown  in  Fig. 

5-11.  The  limit  cycle  behavior  becomes  apparent  in  this 
time  history,  and  as  expected,  the  airplane  remains  in  a 
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Figure 
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Figure  5-10  Throttling  During  Flat  Spin 
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Figure  5-11  Simulation  of  Throttling  in  a  Flat  Spin 


mildly  oscillatory  spin. 

For  an  F-15,  higher  thrust  levels  produce  spin  regions 
with  smaller  stable  branches.  Additionally,  higher  thrust 
levels  produce  spins  at  lower  angles  of  attack  for  a  given 
elevator  deflection.  These  effects  can  aid  in  the  recovery 
from  spins,  however,  used  alone  increasing  thrust  cannot 
bring  the  F-15  out  of  a  spin.  Grafton  (13)  came  to  this 
same  conclusion  in  her  research.  She  determined  that  thrust 
effects  are  generally  small,  but  have  a  generally  favorable 
effect  on  the  number  of  turns  required  to  recover  from  a 
relatively  nonoscil latory  spin. 

Asymmetric  Thrust 

The  F-15  has  the  capability  of  providing  direct  yaw 
moments  by  using  asymmetric  thrust.  In  a  spin,  the  engines 
can  theoretically  be  throttled  to  provide  an  antispin 
moment.  In  real  world  situations,  asymmetric  thrust  usually 
occurs  through  the  inadvertent  flameout  of  one  engine.  Both 
of  these  cases  will  be  looked  at  to  see  how  asymmetric 
thrust  can  lead  to  spins,  and  to  see  how  it  can  be  used  to 
recover  from  them. 

Starting  from  a  low  angle  of  attack  equilibrium  state 
and  a  thrust  of  8300  lbs,  each  engine  was  throttled  and  a  - 
T  bifurcation  diagrams  were  examined.  Figs.  5-12  and  5-13 
show  that  the  equilibrium  branch  continued  to  the  high  angle 
of  attack  regime  through  the  application  of  both  negative 
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and  positive  thrust.  The  negative  thrust  is  unrealistic  but 
does  provide  a  pathway  to  high  a  solutions  in  realistic 
asymmetric  thrust  regions.  The  low  angle  of  attack  branch 
remains  stable  throughout  the  range  of  asymmetries  and 
therefore,  jumps  to  the  spin  region  via  limit  points  should 
not  be  expected  due  to  asymmetric  thrust  alone.  A  large 
perturbation  may  push  the  F-15  away  from  the  stable  branch 
and  cause  an  attraction  to  the  stable  high  a  solutions. 

Fig.  5-12  contains  a  large  region  of  stable  equilibria 
between  a  =  68  degrees  and  80  degrees  and  identifies  a  right 
flat  spin.  This  corresponds  to  the  large  region  of  stable 
right  spin  found  in  the  symmetric  bifurcation  diagrams.  In 
this  region,  T,  >  Tr  and  a  positive  yaw  moment  results.  The 
positive  yaw  moment  drives  the  plane  into  the  right  spin.  A 
left  spin  branch  is  also  present  for  an  asymmetry  of  Tr  >  T, 
but  contains  no  stable  branch.  This  branch  corresponds  to 
the  left  spin  found  in  the  symmetric  bifurcation  diagrams. 
Fig.  5-13  mirrors  Fig.  5-12  but  also  identifies  a  small  sta¬ 
ble  left  spin  branch  at  T,  =  -1000  lbf.  Although  this 
thrust  setting  is  unrealistic  in  the  real  world,  it  shows 
that  an  asymmetry  of  around  5000  lbf  should  be  able  to  pro¬ 
duce  a  stable  left  flat  spin.  The  two  figures  also  identify 
a  small  intermediate  stable  spin  region  between  a  =  50 
degrees  and  58  degrees  that  is  bounded  by  Hopf  points. 
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Left  Engine  =  4150  lbs 
<$.  =  -19.14° 

6r  =  0.0° 


Analysis  of  the  effects  of  asymmetric  thrust  due  to  a 
flameout  were  done  using  elevator  sweeps  for  both  right  and 
left  engine  flameouts  from  an  initial  thrust  of  8300  lbf. 
Figs.  5-14  and  5-15  are  the  a-6t  bifurcation  diagrams  of  a 
right  engine  flameout  and  a  left  engine  flameout.  The  full 
continuation  diagrams  can  be  found  in  Appendix  C.  At  low 
angles  of  attack,  the  onset  of  wing  rock  is  actually  delayed 
by  the  engine  flameout.  Thus,  the  likelihood  of  the  F-15 
jumping  to  the  spin  branch  from  wing  rock  is  not  increased 
by  an  engine  flameout.  However,  Fig  5-14  shows  a  large 
stable  spin  region  extending  from  6.  =  -15  degrees  to  6.  =  13 
degrees.  The  branch  loses  stability  through  Hopf  bifurca¬ 
tions  at  both  ends.  The  positive  yawing  moment  due  to  the 
thrust  asymmetry  creates  a  much  larger  stable  right  spin 
branch  at  higher  values  of  6.  than  the  symmetric  a  -  6, 
bifurcation  diagram.  Additionally,  the  small  left  spin 
branch  virtually  vanishes  from  the  diagram.  Fig.  5-15 
should  produce  an  opposite  effect.  The  stable  right  spin 
branch  stretches  from  6,  =  -26.5  degrees  to  6.  =  -22  degrees 
which  is  smaller  and  at  lower  values  of  6,  than  the  symmet¬ 
ric  diagrams.  The  left  spin  branch  also  becomes  more  preva¬ 
lent.  Thus,  engine  flameouts  can  be  either  helpful  or 
harmful,  depending  on  the  direction  of  spin  and  the  engine 
that  flames  out. 

The  yawing  moment  provided  by  a  thrust  asymmetry  can  be 
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Figure  5-14  Elevator  Sweep  with  Right  Engine  =  0  Lbs 
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used  to  recover  from  spins.  Looking  back  at  Fig.  5-12,  the 
equilibrium  branch  turns  back  when  the  right  engine  thrust 
is  greater  than  15000  pounds.  Any  setting  greater  than 
15000  lbs.  ’’should"  bring  the  F-15  right  back  to -the  stable 
low  angle  of  attack  equilibrium.  A  simulation  was  run  to 
see  the  time  history  of  two  imposed  asymmetries.  The  time 
history  shows  that  the  plane  enters  a  unstable  limit  cycle 
for  a  thrust  setting  of  12000  lbf.  A  full  right  engine 
thrust  setting  of  25,000  lbf  was  then  tried  in  the  simu¬ 
lation.  Fig.  5-16  shows  that  using  the  higher  thrust  asym¬ 
metry,  the  F-15  pushes  through  the  limit  cycle  and  recovers 
to  the  low  angle  of  attack  stable  equilibrium. 
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This  recovery  is  in  a  timeframe  that  could  realistically 
recover  the  P-15  from  the  spin  (13  seconds). 

Asymmetric  thrust  settings  can  aid  in  recovering  from 
spins,  however,  the  same  asymmetries  that  can  aid  the  recov¬ 
ery,  can  make  recovery  more  difficult  if  the  wrong  engine 
flames  out  in  the  spin. 

Thrust  Vectoring 

The  capability  to  vector  thrust  on  an  F-15  is  currently 
available  only  on  the  F-15  STOL  demonstrator.  This  aircraft 
has  modified  nozzles  to  provide  pitch  vectoring,  but  no  yaw 
vectoring  capability  exists.  Pitch  and  yaw  vectoring  will 
be  investigated  to  see  how  they  can  influence  spin  recover¬ 
ies  if  a  baseline  F-15  were  fitted  with  pitch  or  yaw  vector¬ 
ing  nozzles. 

Pitch  Vectoring 

Pitch  vectoring  can  be  thought  of  as  an  automatic  eleva¬ 
tor,  providing  the  pilot  the  capability  to  control  the  pitch 
rate  with  the  engine.  The  pitch  moment  created  by  vectoring 
the  engine  exhaust  is  dependent  on  the  engine  thrust  and  the 
angle  of  vectoring.  Additionally,  vectoring  the  thrust 
decreases  the  force  in  the  x  direction  and  changes  the 
forces  in  the  z  direction.  The  effects  of  pitch  vectoring 
in  aiding  the  recovery  from  spins  are  expected  to  be  small, 
since  the  moment  produced  will  not  oppose  the  large  yaw 
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angular  momentum.  Some  coupling  effects  will  possibly  be 
helpful  in  aiding  the  recovery  and  the  nose  down  moment 
should  lower  the  angle  of  attack  and  make  the  aerodynamic 
controls  more  effective.  Additionally,  a  nose  up  moment 
could  produce  a  deeper  flat  spin. 

An  a  -  6pv  bifurcation  diagram  was  created  for  each  of 
two  different  thrust  levels  (8300  lbs  and  29200  lbs)  start¬ 
ing  from  6,  =  -19.14  degrees  in  the  stable  right  flat  spin 
region.  Figs.  5-17  and  5-18  show  that  a  nose  down  pitching 
moment  (6pv  <  0)  will  cause  the  F-15  to  remain  at  approxi¬ 
mately  the  same  angle  of  attack,  while  a  nose  up  pitching 
moment  (6pv  >  0)  will  bring  about  higher  angles  of  attack. 

The  nose  down  moment  will  also  bring  about  a  loss  of  stabil¬ 
ity  through  a  Hopf  bifurcation.  The  higher  thrust  level 
will  cause  the  limit  cycles  to  appear  at  a  much  lower  pitch 
vector  angle.  The  lower  angle  of  attack  and  loss  of  stabil¬ 
ity  could  possibly  lead  to  recovery.  Additionally,  limit 
cycle  behavior  appears  in  the  nose  up  pitch  vector  for  29200 
lbs.  These  diagrams  were  not  very  helpful  in  trying  to 
understand  the  dynamics  of  the  pitch  vectoring,  so  addi¬ 
tional  diagrams  plotting  p,  q,  and  r  versus  6pv  were  made. 
Figs.  5-19  and  5-20  show  that  coupling  occurs  due  .to  the 
pitch  vectoring.  The  p  -  q  coupling  was  mentioned  in  pre¬ 
vious  discussions  as  the  cause  of  the  6,  anomaly.  Since  p 
is  positive  and  q  is  negative,  the  coupling  term  drives  r 
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into  the  spin  direction.  If  either  term  changed  signs,  the 
coupling  would  provide  an  antispin  term  in  f  and  possibly 
recover  the  F-15  from  the  spin.  Nose  down  pitch  vectoring 
actually  causes  q  to  become  more  negative  and  drives  the 
prospin  yaw  rate  even  higher.  A  nose  up  pit"h  vector  brings 
q  closer  to  zero  and  in  the  case  of  Fig.  5-20  actually  makes 
q  positive.  With  the  positive  q,  an  antispin  coupling 
results  and  the  F-15  should  recover  from  the  spin.  This  is 
not  intuitive  to  a  pilot  since  the  usual  procedure  fcr  spin 
recovery  is  to  release  back  pressure  on  the  stick  which  in 
turn  causes  a  nose  down  pitching  moment.  A  simulation  was 
run  to  see  if  in  fact  this  coupling  takes  place.  Fig.  5-21 
doesn’t  show  that  a  nose  up  moment  will  promote  the  recovery 
from  a  spin  since  the  angle  of  attack  remains  around  90 
degrees.  Although  Fig.  5-21  shows  that  the  F-15  remains  at 
90  degrees,  Fig.  5-22  qualifies  this  by  showing  that  all  yaw 
rotation  stops.  The  F-15  is  in  a  deep  stall  that  can  be 
recovered  by  vectoring  the  nozzle  the  other  direction. 
Applying  nose  down  pitch  vectoring  as  an  initial  command 
actually  increases  the  yaw  rate  and  produces  an  oscillatory 
spin.  Again,  in  the  real  world  an  F-15  wouldn't  have  60 
seconds  to  recover  from  a  flat  spin,  but  these  simulations 
do  show  that  a  nose  up  pitch  vector  can  theoretically  bring 
about  a  recovery.  Pitch  vectoring  is  therefore  useful  in 
recovering  the  F-15  from  spins  but  not  in  an  intuitive  way. 
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Yaw  Vectoring 


Yaw  vectoring  can  also  be  thought  of  as  an  automatic 
control,  providing  the  pilot  control  of  the  yaw  moments  with 
the  engines.  With  yaw  vectoring,  the  capability  exists  to 
directly  oppose  the  angular  momentum  in  a  spin.  Thus,  yaw 
vectoring  is  expected  to  be  the  most  effective  use  of  thrust 
to  recover  from  a  spin.  The  effectiveness  of  yaw  vectoring 
to  recover  from  a  spin  is  directly  proportional  to  the 
thrust  level  and  the  angle  of  deflection.  At  low  thrust 
levels,  the  moment  created  by  vectoring  the  thrust  may  be 
small,  even  with  a  large  vectoring  angle.  At  high  thrust 
levels,  only  a  small  angle  may  be  necessary  to  recover  from 
the  spin. 

The  global  effect  of  yaw  vectoring  was  found  by  starting 
at  a  low  angle  of  attack  equilibrium  point  and  varying  the 
yaw  angle  for  two  thrust  settings.  Figs.  5-23  and  5-24  are 
the  a  -  6yv  bifurcation  diagrams.  The  diagrams  continue  up 
to  the  high  angle  of  attack  regime  where  both  left  and  right 
spin  branches  are  found.  Looking  at  the  low  angle  of  attack 
solutions,  the  equilibrium  solutions  reach  a  turning  point 
at  a  relatively  small  yaw  vector.  The  F-15  may  depart  con¬ 
trolled  flight  at  the  limit  point  and  jump  to  the  spin 
branch  directly  above  it.  Thus,  yaw  vectoring  capabilities 
can  be  dangerous  at  low  angl's  of  attack. 

Looking  at  the  high  a  branches,  the  branch  with  the 
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large  stable  portion  is  the  right  spin  branch.  Since  a 
right  spin  is  the  only  possible  spin  in  this  diagram  for  an 
undeflected  flow,  it  will  be  the  main  topic  of  discussion. 
Both  diagrams  show  that  an  antispin  yaw  vector  (6yv  >0) 
will  drive  the  F-15  out  of  the  stable  spin  and  into  an 
oscillatory  spin  at  the  Hopf  bifurcation.  Further  yaw  vec¬ 
toring  will  bring  the  F-15  to  the  previously  identified 
small  stable  branch  at  a  =  50  degrees.  This  is  surrounded 
by  additional  limit  cycles.  With  enough  yaw  vectoring,  the 
F-15  will  eventually  make  it  through  the  limit  cycles  and 
down  to  the  stable  low  a  equilibrium  branch.  A  prospin  yaw 
vector  will  drive  the  F-15  into  a  deeper,  oscillatory  flat 
spin.  Although  the  diagrams  are  fairly  similar  in  struc¬ 
ture,  the  amount  of  vectoring  necessary  to  drive  the  F-15 
out  of  the  spin  is  much  less  with  higher  thrust  levels.  A 
simulation  was  run  to  see  how  yaw  vectoring  can  bring  the 
F-15  out  of  the  spin.  Fig.  5-25  shows  that  for  a  small  yaw 
vector  (5  degrees),  the  F-15  begins  to  oscillate  around  the 
Hopf  point  discussed  earlier  and  eventually  is  attracted  to 
the  low  a  stable  branch.  With  a  10  degree  yaw  vector,  the 
F-15  is  driven  straight  down  to  the  low  a  stable  branch.  An 
even  larger  yaw  vector  of  20  degrees  brings  the  F-15  down 
to  the  low  a  stable  branch  even  faster,  however,  without 
correction,  the  F-15  will  enter  a  spin  in  the  opposite 
direction.  Looking  back  at  Fig.  5-23  verifies  that  this 
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Thrust  =  8300  lbs 
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Yaw  Thri 

Figure  5-23  Yaw  Vect 
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Yaw  Thrust  Angle 

Figure  5-24  Yaw  Vectoring  T  =  29200  lbs 


would  indeed  happen.  The  low  angle  of  attack  stable  branch 
turns  at  16  degrees  and  the  only  solution  would  be  the  left 
spin.  Fig.  5-26  gives  an  appreciation  as  to  how  quick  the 
rotation  in  the  opposite  direction  begins  if  the  yaw  vector 
is  not  removed  when  the  F-15  approaches  low  angle  of  attack 
flight . 


As  expected,  yaw  vectoring  can  be  very  effective  in 
recovering  the  F-15  from  a  flat  spin.  Burk  (8)  also  came  to 
this  result  in  his  early  research.  However,  yaw  vectoring 
can  lead  to  jumps  to  spins  if  used  at  low  alpha.  Addition¬ 
ally,  care  must  be  taken  to  remove  yaw  vectoring  as  the  F-15 
recovers  or  an  inadvertent  spin  in  the  opposite  direction 
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may  be  entered. 


TIME  (sec) 

Figure  5-26  r  vs  Time  in  Yaw  Vectoring  Simulation 
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Conclusions 


VI . 

Bifurcation  analysis  is  a  powerful  tool  in  the  analysis 
of  nonlinear  aircraft  behavior.  If  it  is  used  in  conjunc¬ 
tion  with  a  realistic  model  of  an  aircraft's  inertia  and 
aerodynamic  coefficients  and  some  trial  simulations,  it  can 
provide  a  qualitative  idea  of  the  nonlinear  behavior  of  the 
aircraft.  From  a  single  starting  point,  it  can  be  used  to 
map  out  an  entire  spectrum  of  possible  aircraft  motions. 

This  study  took  a  fairly  realistic  model  of  an  F-15,  modi¬ 
fied  it  with  nozzles  that  have  never  been  used  on  an  opera¬ 
tional  F-15,  and  mapped  out  how  they  could  be  used  to  help 
recover  from  a  flat  spin.  Additionally,  spin 
characteristics  due  to  thrust  and  thrust  asymmetries  were 
identified.  The  following  conclusions  were  formed  from  this 
research: 

1.  The  F-15  model  designed  by  Baumann  is  a  fairly  real¬ 
istic  model  of  the  F-lE's  actual  behavior.  However,  the 
effectiveness  of  the  rudder  at  high  angles  of  attack  is 
overestimated  in  the  model .  Recovery  from  the  principal 
spin  region  identified  in  the  research  does  not  follow  the 
recommended  course  of  action  of  applying  ailerons  in  the 
direction  of  the  spin.  The  model  therefore  has  some  short¬ 
comings  that  create  unrealistic  aircraft  behavior. 

2.  Thrust  affects  the  spin  characteristics  of  the  F-15 
by  changing  the  size  and  location  of  the  stable  spin 
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branches.  Additionally,  higher  levels  of  thrust  make  thrust 
vectoring  more  effective.  Full  thrust  alone  cannot  recover 
the  F-15  from  a  spin,  but  it  can  aid  other  methods  of  recov¬ 
ering  . 

3.  Asymmetric  thrust  can  be  used  to  recover  the  F-15 
from  a  flat  spin.  The  same  moments  chat  make  thrust  asymme¬ 
tries  useful  can  also  make  spins  more  difficult  to  recover 
from  if  the  wrong  engine  flames  out.  Additionally,  thrust 
asymmetries  do  not  lead  to  jumps  to  spins  and  actually  may 
delay  the  onset  of  wing  rock. 

4.  Pitch  vectoring  can  be  used  in  a  nonintuitive  manner 
to  bring  about  the  recovery  from  a  flat  spin.  Most  pilots 
would  argue  that  a  nose  up  moment  would  probably  deepen  a 
spin,  however  this  research  shows  that  just  the  opposite 
occurs.  Application  of  a  nose  up  moment  reduces  the  yaw 
rate  and  results  in  a  deep  stall  which  can  be  more  easily 
recovered  from. 

5.  Yaw  vectoring  shows  to  be  the  most  promising  method 
of  spin  recovery.  However  care  must  be  taken  to  remove  the 
yaw  vector  upon  recovery  or  an  inadvertent  spin  in  the  oppo¬ 
site  direction  can  occur.  The  use  of  yaw  vectoring  at  low  a 
can  lead  to  jumps  to  spins  due  to  limit  point  behavior  and 
therefore  should  be  avoided. 

Pecommendations 
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The  following  are  recommendations  for  future  research 
that  were  identified  while  pursuing  this  research: 

1.  Use  rotary  balance  data  in  the  model.  Many  previous 
studies  of  spin  behavior  used  rotary  balance  data  in  their 
models.  F-15  rotary  balance  data  is  available  but  was  not 
used  in  this  study  (3,4).  This  data  is  nonlinear  in  angular 
velocity  and  more  correctly  characterizes  the  actual  aerody¬ 
namics  in  a  spin.  Discrepancies  between  flight  test  data 
and  predicted  behavior  identified  in  this  research  may  be 
explained  away  using  the  rotary  balance  data. 

2.  Plug  the  thrust  vectoring  angles  into  a  control  sys¬ 
tem  to  see  how  it  can  be  used  to  minimize  wing  rock.  This 
offers  an  ideal  method  to  help  keep  control  of  the  P-15  at 
high  angles  of  attack  by  not  having  to  rely  on  aerodynamic 
surfaces . 

3-  The  F-15  used  in  this  study  was  loaded  symmetri¬ 
cally.  In  the  real  world,  the  F-15  is  often  flown  with  an 
asymmetric  load  of  fuel  or  the  weapons.  This  is  a  potential 
cause  of  many  departures  in  the  F-15  and  bifurcation  theory 
can  be  used  to  better  characterize  the  mechanisms  of  these 
departures.  Guicheteau  (16:7)  has  applied  this  in  his  study 
of  the  Alpha  Jet. 

4.  Include  the  contribution  of  gyroscopic  torques  due 
to  the  engines.  This  analysis  assumes  that  they  are  zero. 

In  reality,  they  will  influence  the  spin  behavior  of  the 
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F-15.  Guicheteau  (16:7)  also  applies  this  to  the  Alpha  Jet. 

5.  The  model  has  an  inaccuracy  identified  (Cn8r)  and  may 
contain  more  that  were  not  identified  in  this  research. 
Therefore,  a  thorough  refinement  of  the  model’s  curve  fits 
could  make  it  more  realistic  in  the  high  a  regime. 

6.  Apply  bifurcation  theory  to  one  of  the  new  designs 
(ATF,  B-2,  C-17)  using  wind  tunnel  data  to  help  develop 
their  flight  test  program  or  verify  flight  test  results. 
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Appendix  A:  F-15B  Weight  and  Balance  Data 


The  physical  dimensions  and  weight  and  balance  data  for 
the  F-15B  is  listed  in  Table  V.  The  data  for  the  F-15  was 
obtained  from  Beck  (7)  and  (23). 

Table  II.  Physical  Characteristics  of  the  F-15B 


Wing 

Area  (Theoretical) 

Aspect  Ratio 
Airfoil 
Root 
Xw  155 
Tip 
Span 

Taper  Ratio 

Root  Chord  (Theoretical) 
Tip  Chord 

Mean  Aerodynamic  Chord 

Leading  Edge  Sweep  Angle 

25%  Chord  Sweep  Angle 

Dihedral 

Incidence 

Twist  at  Tip 

Aileron  Area 

Flap  Area 

Speed  Brake  -  Area 

Control  Surface  Movement 
Aileron 
Speedbrake 
Flap 

Horizontal  Tail 
Rudder 

Vertical  Tail 

Area  (Theoretical  Each) 
Rudder  Area  (Each) 

Span 

Aspect  Ratio 
Root  Chord 
Tip  Chord 
Airfoil  -  Root 


608  sq  ft 
3.01 

NACA64006 . 6 

NACA64A( x)04 . 6  (a  =  0.8  Mod) 
NACA64A203  (a=0.P  Mod) 

42.8  ft 
0.25 

273.3  in 

68.3  in 

191.3  in 
45  degrees 

38.6  degrees 
-1  degrees 
None 

None 

26.5  sq  ft 

35.8  sq  ft 

31.5  sq  ft 


+/-  20  degrees 
45  degrees  up 
30  degrees  down 

29  degrees  down,  15  degrees  up 
+/-  30  degrees 


62.6  sq  ft 
10.0  sq  ft 

10.3  sq  ft 
1.70 

115.0  in 

30.6  in 
NACA0005-64 


69 


-  Tip 

NACA0003 . 5-64 

Taper  Ratio 

0.27 

Leading  Edge  Sweep  Angle 

36.6  degrees 

25%  Chord  Sweep  Angle 

29.7  degrees 

Mean  Aerodynamic  Chord 

81.0  in 

Cant 

2  degrees  out 

Length  (.25cwto  .25ch) 

241.0  in 

Wetted  Area 

Fuselage 

1405  sq  ft 

Nozzles 

53  sq  ft 

Horizontal  Tail 

216  sq  ft 

Vertical  Tail 

257  sq  ft 

Wing 

698  sq  ft 

Total  Area 

2629  sq  ft 

Engine  Data  (each) 

Non-Afterburning  Thrust 

14,871  lb 

Afterburning  Thrust 

23,810  lb 

Y  Direction  C.G.  Offset 

+/-  25.5  in 

Z  Direction  C.G.  Offset 

0.25  in 

Nozzle  Pivot  C.G.  Offset 

-20.219  ft 

Miscellaneous  Data 

Aircraft  Length 

63.8  ft 

Aircraft  Height 

18.6  ft 

Aircraft  Volume 

1996  cu  ft 

Aircraft  Gross  Weight 

37000  lbs 

C.G.  Station  X  Direction 

557.173 

Y  direction 

0.0 

Z  Direction 

116.173 

Inertia  Data 

I* 

25480  slug-ft2 

Iy 

166620  slug-ft 

U 

186930  slug-ft 

lx. 

-1000  slug-ft2 

The  inertia  values  are  for  a  basic  F-15  with  4  AIM- 
missiles,  ammo,  50%  fuel  and  gear  up. 
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Appendix  B:  Driver  Program 


C  CAPTAIN  ROBERT  J.  MCDONNELL  AFIT  GAE-90D 

C  MASTERS  THESIS 

C 

C  THIS  COMPUTER  PROGRAM  SOLVES  THE  NONLINEAR  DIFFERENTIAL 
C  EQUATIONS  OF  MOTION  FOR  THE  F-15B  AIRCRAFT.  IT  IS  USED 
C  AS  AN  ANALYTICAL  TOOL  IN  THE  SEARCH  OF  HIGH  ANGLE  OF  ATTACK 

C  PHENOMENA  (I.  E.  FLAT  SPINS).  THE  PROGRAM  IS  CAPABLE  OF 

C  VARYING  ELEVATOR,  AILERON,  AND  RUDDER  DEFLECTIONS,  ENGINE 
C  THRUST  VECTOR  (PITCH  AND  YAW) ,  PORT  AND  STARBOARD  ENGINE 
C  THRUST,  AND  TOTAL  THRUST. 

C 

C  LAST  EDITED  ON  24  OCT  1990 

C 

IMPLICIT  DOUBLE  PRECISION (A-H,0-Z) 

DIMENSION  W( 300000),  IW(1000) 

OPEN(UNIT=3 ,FILE= ' fort . 3 ' ) 

OPEN(UNIT=4 , FILE= ' f ort . 4 ' ) 

OPEN(UNIT=7 , FILE= ' fort . 7 ' ) 

OPEN(UNIT=8 , FILE= ' fort .  8 ' ) 

OPEN (UNIT= 9 , FILE= ' fort . 9 ' ) 

OPEN(UNIT=10 ,FILE= ' fort . 10 ' ) 

OPEN (UNIT=12 , FILE= ' cs 1 ) 

OPEN (UNIT=13 , FILE= ’ cts ’ ) 

REWIND  7 
REWIND  8 
REWIND  9 
REWIND  10 
REWIND  3 
REWIND  4 
REWIND  12 
REWIND  13 

CALL  AUTO  -  CONTINUATION  &  BIFURCATION  LOCATION  SUBROUTINE 

CALL  AUTO(W,IW) 

STOP 
END 

SUBROUTINE  FUNC(NDIM,NPAR,U,ICP,PAR, IJAC,F,DFDU,DFDP) 


IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

COMMON  /KS/  K1,K5,K7,K8,K9,K10,K12,K13,K14,K15,K16,K17 
COMMON  /ACDATA/  SWING, OWING, SREF,RHO,RMASS 
DOUBLE  PRECISION  K1,K5,K7,K8,K9,K10,K12,K13,K14,K15,K16,K17 
CCM-iON  /SEIZE/  CX,CY,CZ ,CIM,CM4,CNM 
CCM40N  /SEIZET/  CXT,CYT,CZT,CI2fr,C»lT,CNMT 
C 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


DIMENSION  DFDU(NDIM,NDIM)  ,DFDP(NDIM,NPAR)  ,DELF1(8) , 

+  DE1iF2(8)  ,U(NDIM)  ,PAR(10)  ,F(NDIM)  ,DX(8) 

INITIALIZE  SOME  CONSTANTS  THAT  ARE  PASSED  THROUGH 
THE  COMMON  BLOCK  ACDATA 

DATA  IS  FROM  MCAIR  REPORT#  A4172  AND  AFFTC-TR-75-32 
F-15A  APPRQACH-TO- STALL/ STALL/POST- STALL  EVALUATION 

SWING  -  A/C  WINGSPAN,  FT 

OWING  -  A/C  MEAN  AERODYNAMIC  CHORD,  FT 

SREF  -  A/C  WING  REFERENCE  AREA,  SQ  FT 

RHO  -  AIR  DENSITY  AT  20000  FT  ALTITUDE,  SLUG/FT3 

RMASS  -  A/C  MASS,  SLUGS 

BWING=42.8 

CWING=15.94 

SREF=608. 

RHO=. 0012 673 
RMASS= 37000 . /32 . 174 

DETERMINE  CONSTA!JTS  K1  THROUGH  K17.  SOME  ARE  MADE  COMMON  AND 
PASSED  TO  SUBROUTINE  FUNX  AND  USED  IN  THE  EQUATIONS 
OF  MOTION  THERE 

INERTIAS  HAVE  UNITS  OF  SLUG-FT2 
K1  HAS  UNITS  OF  1/FT 

K6,  K8,  Kll,  K14,  AND  K17  HAVE  UNITS  OF  1/FT2 

IX=  25480. OdO 

IY=  166620. OdO 

IZ=  186930. OdO 

IXZ=  -1000. OdO 

K1=0 . 5dO*RHO*SREF/RMASS 

K2=(IZ-IY)/IX 

K3=IXZ*IXZ/(IX*IZ) 

K4=(IY-IX)/IZ 

K5=IXZ/IX 

K6=0 . 5dO*RHO*BWING*SREF/IX 
K7=IXZ/IZ 

K8=0 . 5dO*RHO*SREF* OWING/ IY 

K9=(IZ-IX)/IY 

K10=IXZ/IY 

K11=0 . 5dO*RHO*SREF*BWING/IZ 
K12= (K2+K3 ) / ( 1 . 0d0-K3 ) 

K13=(1.0d0-K4)*K5/(1.0d0-K3) 

K14=K6/ ( 1 . OdO -K3 ) 

K15= (K3-K4 ) / ( 1 . 0d0-K3 ) 

K16=(1.0d0+K2)*K7/(1.0d0-K3) 

K17=Kll/(1.0d0-K3) 
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K1  =  3 . 350088890D-04 
K5  =-3 . 924646781D-02 
K7  =-5 . 349596105D-03 
K8  =  3.685650971D-05 
K9  =  .96897131196 
K10  =-6 . 001680471D-03 
K12  =  .79747314581 
K13  =-9.615755341D-03 
K14  =  6.472745847D-04 
K15  =-.754990553922 
K16  =  K13 

K17  =  8 . 822851558D-05 
C 

C  FIND  THE  VALUES  OF  F(l)  THROUGH  F(NDIM) .  SUBROUTINES 
C  COEFF  AND  FUNX  ARE  CALLED  CWCE. 

C 

CALL  COEFF(U,PAR,NDIM,ICP) 

C 

CALL  roNX(NDIM,UfF) 

C 

IF(IJAC.EQ.O)  RETURN 
C 

C  SET  THE  VALUES  OF  DX 

C  MODIFIED  TO  SCALE  DX  ACCORDING  TO  VARIABLE 
C  13  JUN  88 
C 

C  DX0=1.0D-9 

DX( 1)=DX0*50 . OdO 
DX( 2 ) =DX0*10 . OdO 
DX( 3 ) =DX0*0 . 5d0 
DX( 4 ) =DX0*0 . 25d0 
DX( 5 ) =DX0*0 . 5d0 
DX ( 6 ) =DX0*50 . OdO 
DX(7 )=DX0*50 . OdO 
DX( 8 ) =DX0*0 . 5d0 
C 

C  NEXT  THE  PARTIAL  OF  F  W.R.T.  A  GIVEN  PARAMETER  ARE  FINITE 
C  DIFFERENCED 
C 

PTEMP=PAR(ICP) 

PAR(ICP)=PTEMP+DX(1) 

CALL  COEFF(U,PAR,NDIM, ICP) 

CALL  FUNX(NDIM,U,DELF1) 

C 

PAR( ICP)=PTEMP-DX(1) 

CALL  COEFF(U,PAR,NDIM,ICP) 

CALL  FUNX(NDIM,U,DELF2) 

C 

DO  13  I=1,NDIM 
C 

DFDP( I , ICP)= (DELF1 ( I ) -DELF2 ( I ) ) / ( 2 . OdO*DX( 1 ) ) 

C 

13  CONTINUE 
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PAR(ICP)=PTEMP 


THE  NEXT  DO  LOOP  CALCULATES  THE  PARTIAL  DERIVATIVE  OF  F  W.R.T. 
TO  U  USING  FINITE  DIFFERENCES. 

SET  U(J)  EQUAL  TO  U+DU,  THEN  CALL  COEFF  WITH  THIS  UPDATED 
STATE  VECTOR.  THIS  IS  DONE  SIMILARLY  WITH  U-DU 

DO  20  J=1,NDIM 

UTEMP=U(J) 

U(J)=UTEMP+DX(J) 

CALL  COEFF(U,PAR(NDIM,  ICP) 

CALL  EUNX(NDIM,U,D£LF1) 

U(J)=UTEMP-DX(J) 

CALL  OOEFF(U, PAR, NDIM, ICP) 

CALL  EUNX(NDIM,U,DELF2) 

DO  16  I=1,NDIM 

DFDU( I , J)= (DELF1( I ) -DELF2( I ) )/ (2 . 0d0*DX( J) ) 

16  CONTINUE 

U(J)=UTEMP 

20  CONTINUE 
RETURN 
END 

SUBROUTINE  FUNX(NDIM,U,F) 


SUBROUTINE  FUNX  EVALUATES  THE  NDIM  EQUATIONS  GIVEN  THE 
STATE  VECTOR  U. 

NDIM-  THE  DIMENSION  OF  THE  PROBLEM 
U  -  THE  VECTOR  OF  STATES  ALPHA,  BETA,  ...  (INPUT) 

F  -  THE  VECTOR  RESULT  OF  FUNCTION  EVALUATIONS  (OUTPUT) 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

COMMON  /SEIZE/  CX,CY,CZ,CLM,C**4,CNM 
COMMON  /SEIZET/  CXT,CYT,CZT,CLMr,CJ<Wr,CNMr 
COMMON  /KS/  K1,K5,K7 ,K8 ,K9,K10,K12 ,K13,K14,K15,K16,K17 
DOUBLE  PRECISION  K1,K5,K7,K8,K9,K10,K12,K13,K14,K15,K16,K17 
DIMENSION  U(NDIM)  ,F(NDIM) 

SET  TRIGONOMETRIC  RELATIONSHIPS  OF  THE  STATES  ALPHA,  BETA, 
THETA,  AND  PHI  AND  THEN  SET  P,  Q,  R,  AND  VTRFPS 

IWRITE=1 
C 

DEGRAD=57 . 29577951D0 
C 


74 


ooooo  OOOOO  ooooo  ooooooooooo  o  o  o 


CA=C0S(U(1)/DBSRAD) 

SA=SIN(U(1)/DEGRAD) 

CB=COS(U(  2) /DEGRAD) 

SB=SIN(U(2)/DEGRAD) 

cthe=oos(u(6)/dbgrad) 

STHE=SIN(U(6)/DEGRAD) 

CFHI=OOS(U(  7) /DEGRAD) 

SPHI=SIN(U(7)/DEGRAD) 

P=U(3) 

0*0(4) 

R=U(5) 

VTRFPS=1000.0d0*U(8) 

SET  THE  GRAVITATIONAL  CONSTANT,  FT/ SEC 
G=32.1740d0 

THE  FOLLOWING  SYSTEM  OF  NONLINEAR  DIFFERENTIAL  EQUATIONS 
GOVERN  AIRCRAFT  MOTION 

UPDATED  FOR  PROPER  DEGREE-RADIAN  UNITS  AND  PROPERLY 
SCALED  VELOCITY  EQUATION:  7  JUN  88 

******  alpha  dot  ****** 

F ( 1 ) =ALFHA~DOT 

1  F(  1 )  =Q+  (  -  (K1*VTRFPS*CX-G*STHE/VTRFPS+R*SB)  *SA+  (K1*VTRFPS 

+  *CZ+ ( G*CTHE*CFHI / VTRFPS ) -P*SB ) *CA) /CB 

F(1)=F(1)*DEGRAD 

******  BE7TA  DOT  ****** 

F(2)=BETA-D0T 

2  F(2)=-( ( K1*VTRFPS*CX-G*STHE/ VTRFPS ) *SB+R ) *CA+ (  K1*VTRFPS*CY 

+  +G*CIHE*SPHI /VTRFPS )  *CB-  (  (K1*VTRFPS*CZ+G*CTHE*CPHI/VTRFPS) 

+  *SB-P)*SA 

F(2)=F(2)*DEGRAD 

******  p  ixyr  ****** 

F(3)=P-D0T 

3  F(  3  )=-K12*Q*R+K13*P*Q+K14*  (OM+K7*CNM)  *VTRFPS*VTRFPS 
******  q  COT  ****** 

F(4)=Q-D0T 

4  F^4)=K8*VTOFPS*VTRFPS*Ofl+K9*P*R+K10*(R*R-P*P) 

C 
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******  p  J3QT  ****** 

F(5)=R-D0T 

5  F(  5  )  =K15*P*Q-K16*Q*R+K17  *VTRFPS*VTRFPS*  (K5*CU4+CNM) 
******  ihETA  DOT  ****** 

F(6)sTHETA-DCT 

6  F(6)=Q*CFHI-R*SFHI 
F(6)=F(6)*DEGRAD 

******  pjji  ix)T  ****** 

F(7)=PHI-D0T 

7  F(7)=P+Q*(STHE/CniE)*SPHI+R*(STOE/CTHE)*CEHI 
F(7)=F(7)*DBGRAD 

******  v  DOT  ****** 


F ( 8 ) = VTRFPS-DOT  (SCALED  BY  A  FACTOR  OF  1000) 

8  F(  8 )  =U(  8 )  *  ( ( K1*VIRFPS*CX-G*  STHE/ VTRFPS  )  *CA*CB+  (K1*VTRFPS*CY 

+  +G*CTHE*SPHI/VTRFPS)*SB 
+  +  ( Kl*VTRFPS*CZ*3*CTHE*CPHI  /VTRFPS )  *SA*CB ) 

RETURN 

END 

SUBROUTINE  STPNT(NDIM,U,NPAR,ICP,PAR) 


THIS  SUBROUTINE  SETS  THE  VALUES  OF  THE  STATES  AND  PARAMETERS 
AT  THE  START  OF  THE  ANALYSIS.  THE  STATES  AND  CONTROL  SURFACE 
SETTINGS  REPRESENT  AN  EQUILIBRIUM  STATE  OF  THE  AIRCRAFT 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 


IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

DIMENSION  U(NDIM) ,PAR(10) 

U(l)  -  ALPHA,  DEG 
U(2)  -  BETA,  DEG 
U(3)  -  P,  RAD/SEC 
U(4)  -  Q,  RAD/SEC 
U(5)  -  R,  RAD/SEC 
U(6)  -  THETA,  DEG 
U(7)  -  PHI,  DBG 

U(8)  -  TRUE  VELOCITY,  IN  THOUSANDS  OF  FT/SEC 
THE  STARTING  POINT  (VECTOR) 

OPEN(UNIT=15 ,FILE= ’ f ort . 15 ' ) 
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REWIND  (15) 

C 

READ(15,*)  U(l) 

READ(15,*)  U(2) 

READ(15,*)  U(3) 

READ(15,*)  U(4) 

READ(15,*)  U(5) 

READ(15,*)  U(6) 

READ(15,*)  U(7) 

READ(15,*)  VTRFPS 
U(8)=VTRFPS/1000 . OdO 
C 

C  PAR(1)=DELESD 

C  PAR(2)=DRUDD  THE  PARAMETERS,  IN  DEGREES 
C  PAR(3)=DDA 

C  PAR(4)=ENGPA  PORT  ENGINE  THRUST,  POUNDS/ 1000 
C  PAR(5)=ENGSA  STARBORD  ENGINE  THRUST,  POUNDS/IOOO 
C  PAR(6)=TPTAL  PITCH  THRUST  VECTOR,  DEG 

C  PAR(7)=TYTAL  YAW  THRUST  VECTOR,  DEG 

C  PAR(8)=TTHRST  TOTAL  THRUST,  POUNDS/ 1000 

C 

READ(15,*)  PAR(l) 

READ(15,*)  PAR(2) 

READ(15,*)  PAR(3) 

READ(15,*)  PAR(4) 

READ(15,*)  PAR(5) 

READ(15,*)  PAR(6) 

READ(15,*)  PAR(7 ) 

READ(15,*)  PAR(8) 

C 

RETURN 

END 

C 

SUBROUTINE  IN IT 

C  - 

C 

IMPLICIT  DOUBLE  PRECISI0N(A-H,O-Z) 

C 

COMMON  /BLCSS/  NDIM,ITMX,NPAR,ICP,IID,NMX,IPS,IRS 
CttMON  /BLCPS/  NTST,NCOL,IANCH,NMXPS,IAD,NFR,NWTN,ISP,ISWl 
CCM40N  /BLDLS/  DS,D34IN,Da4AX,IADS 
COMMON  /BLLIM/  RL0,RL1,A0,A1,PAR(10) 

COMMON  /BLOPT/  ITNW,MXBF,IPLT,ICP2,ILP 
OQMMDN  /BLEPS/  EPSU , EPSL , EPSS , EPSR 
C 

C  IN  THIS  SUBROUTINE  THE  USER  SHOULD  SET  THOSE  CONSTANTS 
C  THAT  REQUIRE  VALUES  DIFFERENT  FROM  THE  DEFAULT  VALUES 
C  ASSIGNED  IN  THE  LIBRARY  SUBROUTINE  DFINIT.  FOR  A  DESCRIPTION 
C  OF  THESE  CONSTANTS  SEE  THE  DOCUMENTATION  CONTAINED  IN  THE 
C  LIBRARY.  CCM40N  BLOCKS  CORRESPONDING  TO  CONSTANTS  THAT  THE  USER 
C  WANTS  TO  CHANGE  MUST  BE  INSERTED  ABOVE.  THESE  CC*MON  BLOCKS 
C  SHOULD  OF  COURSE  BE  IDENTICAL  TO  THOSE  IN  DFINIT. 

C 
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DSMAX  =  lO.OdO 
D31IN  =0.00000010d0 
EPSU  =  1.0D-07 
EPSL  =  1.0D-07 
EPSS  =  1.0D-05 
EPSR  =  1.0D-07 
I  AD  =1 

ILP  =1 
I7MX  =  40 
ITNW  =  20 
MXBF  =  5 
NDIM  =  8 
NPAR  =  8 
C 

OPEN(UNIT=25 ,FILE= ’ fort . 25 ' ) 

REWIND  (25) 

C 

READ(25,*)  RL0,RL1 
READ(25,*)  A0,A1 
READ(25,*)  DS 
READ(25,*)  NMX 

READ(25,*)  NTST,NCOL,NMXPS,NPR 
READ(25,*)  ISP , IRS , ICP / ICP2 , IPLT , IPS 
READ(25,*)  ISW1 
RETORN 
END 
C 

SUBROUTINE  BOND 


RETORN 

END 

SUBROUTINE  ICND 


RETORN 

END 

SUBROUTINE  COEFF(U, PAR, NDIM,  ICP) 


IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

COMMON  /ACDATA/  BW ING,  OWING,  SREF,RHO,RMASS 
CCM40N  /SEIZE/  CX,CY,CZ,CLM./aM,CNM 
(XM10N  /SEIZET/CXT,CYT,CZT,CLMT,CMfr,ClWr 
DIMENSION  U(NDIM)  ,PAR(10) 

C 

C  TOE  PRIMARY  SOURCE  OF  THESE  COEFFICIENT  EQUATIONS  IS  SUBROUTINE 
C  AR010  FROM  MCAIR  CODE  USED  IN  TOE  F15  BASELINE  SIMULATOR. 

C 

C  MOST  OF  TOE  COEFFICIENTS  USED  IN  TOE  EQUATIONS  WERE  COMPUTED 
C  USING  SAS  WITH  RAW  DATA  FROM  TOE  F15  SIMULATOR  DATA  TABLES. 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


THIS  SUBROUTINE  IS  CALLED  BY  THE  DRIVER  PROGRAM  FOR  THE  AUTO 
SOFTWARE.  IT  MERELY  TAKES  INPUTS  ON  THE  A/C  STATE,  CONTROL 
SURFACE  POSITIONS,  AND  THRUST  VALUES  AND  RETURNS  THE 
APPROPRIATE  AERO  COEFFICIENTS  CX,  CY,  CZ,  CL,  CM,  AND  CM. 

INPUTS  TO  THIS  SUBROUTINE 

AL  -  ANGLE  OF  ATTACK,  DEG 

BETA  -  SIDESLIP  ANGLE,  DEG 

DDA  -  AILERON  DEFLECTION  ANGLE,  DEG 

DELEDD  -  DIFFERENTIAL  TAIL  DEFLECTION  ANGLE,  DEG 

DELESD  -  SYMMETRICAL  TAIL  DEFLECTION  ANGLE,  DEG 

DRUDD  -  RUDDER  DEFLECTION,  POSITIVE  TRAILING  EDGE  LEFT,  DEG 

P  -  ROLL  RATE,  RAD/SEC 

Q  -  PITCH  RATE,  RAD/SEC 

R  -  YAW  RATE,  RAD/SEC 

ENGPA  -  PORT  ENGINE  THRUST,  POUNDS/ 1000 

ENGSA  -  STARBOARD  ENGINE  THRUST,  POUNDS/ 1000 

TYTAL  -  YAW  THRUST  ANGLE,  DEG 

TPTAL  -  PITCH  THRUST  ANGLE,  DEG 

TTHRST  -  TOTAL  THRUST,  POUNDS/IOOO 

VTRFPS  -  TRUE  AIRSPEED,  FT/ SEC 

INTERMEDIATE  VARIABLES  USED  IN  THIS  SUBROUTINE 

ABET  -  ABSOLUTE  VALUE  OF  BETA,  DEG 

ARUD  -  ABSOLUTE  VALUE  OF  RUDDER  DEFLECTION,  DEG 

SWING  -  WING  SPAN,  FEET 

CA  -  COSINE  RAL  (RAL  IN  RADIANS) 

CD  -  COEFFICIENT  OF  DRAG 

CL  -  BASIC  LIFT  COEFFICIENT 

OWING  -  MEAN  AERODYNAMIC  CHORD,  FEET 

DAHD  -  DIFFERENTIAL  ELEVATOR  DEFLECTION,  DEG 

DAHLD  -  LEFT  AILERON  DEFLECTION,  DEG 

DAHRD  -  RIGHT  AILERON  DELFECTION,  DEG 

DELEDR  -  DIFFERENTIAL  TAIL  DEFLECTION  ANGLE,  RAD 

DELESR  -  SYMMETRIC  TAIL  DEFLECTION  ANGLE,  RAD 

ENGP  -  PORT  ENGINE  THRUST,  POUNDS 

ENGS  -  STARBOARD  ENGINE  THRUST,  POUNDS 

PTAL  -  PITCH  THRUST  VECTOR,  RAD 

QBARS  -  DYNAMIC  PRESSURE  TIMES  WING  REFERENCE  AREA,  LBF 

RABET  -  ABSOLUTE  VALUE  OF  BETA,  RADIANS 

RAL  -  ABSOLUTE  VALUE  OF  ALPHA,  RADIANS 

RARUD  -  ABSOLUTE  VALUE  OF  RUDDER,  RADIANS 

SA  -  SINE  RAL  (RAL  IN  RADIANS) 

YTAL  -  YAW  THRUST  VECTOR,  RAD 

OUTPUTS  FRCM  THIS  SUBROUTINE 

CX  -  BASIC  AXIAL  FORCE  COEFFICIENT,  BODY  AXIS,  +  FORWARD 
CY  -  BASIC  SIDE  FORCE  COEFFICIENT,  BODY  AXIS,  +  RIGHT 

CZ  -  BASIC  NORMAL  FORCE  COEFFICIENT,  BODY  AXIS,  +  DOWN 
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C  CLM  -  BASIC  ROLLING  MCMENT  COEFFICIENT,  BODY  AXIS,  +  R  WING 
C  DOWN 

C  CMM  -  BASIC  PITCHING  MCMENT  COEFFICIENT,  BODY  AXIS,  +  NOSE  UP 

C  CNM  -  BASIC  YAWING  MCMENT  COEFFICIENT,  BODY  AXIS,  +  NOSE  RIGHT 

C 

C  ANGLES  USED  IN  CALCULATING  CL,  CLLDB,  . . . ,  ARE  IN  RADIANS.  THIS 
C  IS  BECAUSE  RADIANS  WERE  USED  IN  TOE  CURVE  FITTING  PROGRAM  TO 

C  OBTAIN  TOE  COEFFICIENTS  OF  TOE  ALPHA,  BETA,  . . . ,  TERMS  IN  TOE 

C  FOLLOWING  EQUATIONS. 

C 

C  MCMENT  REFERENCE  CENTER  WAS  SET  IN  AROIO  PROGRAM  AS: 

C 

C  DATA  CMOGR  /.2565/,  CNCGR  /.2565/ 

C 

C  THE  AERO  STABILITY  DATA  WAS  TAKEN  REFERENCED  TO  TOESE  CG 
C  LOCATIONS.  TOE  MOMENTS  OF  INERTIA  AND  OTHER  AIRCRAFT  DATA 

C  ARE  FOR  A  CLEAN  CONFIGURATION  TEST  AIRCRAFT  WITH  A  CG  AT 

C  THE  SAME  CG.  AS  A  RESULT,  THERE  IS  NO  ’CG  OFFSET'  TO  BE 
C  COMPUTED. 

C 

IWRITE=0 

C 

AL=U(1) 

BETA=U(2) 

P=U(3) 

Q=U(4) 

R=U(5) 

TOETA=U(6) 

FHI=U(7) 

VTRFPS=U(8)*1000. 

C 

DELESD=PAR(1) 

DRUDD=PAR(2) 

DDA=PAR(3) 

ENGPA=PAR(4) 

ENGSA=PAR(5) 

TPTAL=PAR(6) 

TYTAL=PAR(7) 

TTHRST=PAR(8) 

C 

DEX3RAD=57. 29577951 
DELESR=DELESD/ DEGRAD 
YTAL=TYTAL/ DEGRAD 
PTAL=TPTAL/DEGRAD 
C 

C  IF  BLOCK  TO  CHANGE  TOTAL  THRUST 
C 

IF(ICP.EQ.8)TOEN 

DIFT=PAR(4)-PAR(5) 

TOALF=TTORST/ 2 . OdO 
ENGPA=THALf +DIFT/2 . OdO 
ENGSA=TOALF-DIFT/2 .  OdO 
ENDIF 
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c 


ENGP=ENGPA*1000.0 

ENGS=E2iGSA*1000.0 


C 

QBARS=0 .  MO*RHO*VTRFPS*VTRFPS*SREF 
C02V=CWING/ ( 2 . OdO*VTRFPS ) 

B02V=BMING/ ( 2 . OdO*VTRFPS ) 

QSB=BWING*QBARS 
ARUD= ABS ( DRUDD ) 

RARUD=ARUD/ DEGRAD 
RAL=AL/ DEGRAD 
ABET=ABS(BETA) 

RABET= ABET/DEGRAD 
C 

(^*********************jlt*********<f************************************* 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 


NEW  SECTION  OF  CODE  -  1)  ALL  THE  AERODYNAMIC  COEFFICIENTS  IN 

THIS  VERSION  OF  TOE  DRIVER  PROGRAM 
ARE  TAKEN  DIRECTLY  ERCM  THE  1988 
F15  AEROBASE  (0.6  MACH,  20000  FEET) 

2)  THIS  SECTION  SUMMARIZES  THE 

AERODYNAMIC  COEFFICIENTS  AS  TO  WHAT 
THEY  ARE  AND  HOW  THEY  ARE  USED. 

TOE  FIRST  ACCRONYM  IS  TOE  JOVIAL  NAME 
OF  TOE  AERODYNAMIC  COEFFICIENT  (CFX1, 
ETC) ,  TOE  SECOND  ACCRONYM  IS  THE 
F15  AEROBASE  CODE  OR  CTAB  NAME 
(ATAB15,  ETC).  A  BRIEF  DEFINITION 
OF  THE  AERODYNAMIC  COEFFICIENT  IS  ALSO 
PROVIDED. 


C 

C 

C 

C 

C 

C 

c 

c 

c 

c 


3)  THERE  IS  ALSO  A  SECTION  THAT  PROVIDES 
A  TABLE  OF  CONVERSIONS  BETWEEN  WHAT 
THE  VARIABLE  IS  CAT  .LED  IN  THE  ORIGINAL 
SECTION  OF  THIS  PROGRAM 
AND  ITS  NAME  IN  TOE  1988  F15  AEROBASE. 
FOR  THE  SAKE  OF  CONTINUITY  THE 
ORIGINAL  PROGRAM  NAME  IS  USED  AND 
THE  1988  F15  AEROBASE  NAME 
IS  PROVIDED  AS  BOOK  KEEPING 
INFORMATION. 


C*********  CFyi  *****************************************************c 


C 

C  CFX  =  FORCE  IN  STABILITY  AXIS  X  DIRECTION  (CD  IN  BODY  AXIS) 

C  (FUNCTION  OF  CL  OR  CFZ1) 

C  CFX  =  CFX1  +  CXRB  +  STORE  INCREMENTS  +  CXDSPD  +  DCXLG  +  DCD 
C 

C  CFX1  =  ATAB15  =  PERFORMANCE  DRAG  COEFFICIENT  -  CD 
C  CXRB  =  ATAB22  =  DELTA  CD  DUE  TO  OG  (=0.0) 

C  CXDSPD  =  ATAB27  =  DELTA  CD  DUE  TO  SPEEDBRAKE  (NORMALLY  =  0.0436) 
C  SET  TO  0  SINCE  THIS  STUDY  IS  CONCERNED 
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c 

c 

c 

c 

c 

c 

c 

c 


WITH  HIGH  ANGLES 

OF  ATTACK  PHENOMENON  (>40  DEGREES)  AND  BECAUSE 
THE  SPEEDBRAKE  WILL  NOT  DEPLOY  AT  ANGLES  OF 
ATTACK  GREATER  THAN  15  DEGREES. 

DCXLG  =  ATAB19  =  DELTA  CD  DUE  TO  REYNOLD'S  NUMBER  (=-0.0005) 

DCD  =  BTAB03  =  DELTA  CD  DUE  TO  2-PLACE  CANOPY  (F15B)  (=0.0005) 

*******  NOTE  THAT  DCXLG  AND  DCD  CANCEL  EACH  OTHER  ******* 


C***********  Qpy  *************************************************q 


C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

n 

V-' 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


CFY  =  FORCE  IN  BODY  AXIS  Y  DIRECTION 

CFY  =  CFY1*EPA02  +  CYDAD*DAILD  +  [CYDRD*DRUDD*DRFLX5]*EPA43 
+[CYDTD^DTFLX5  +  DTFLX6]*DTALD  +  CFYP*PB  +  CFYR*RB 
+CYRB  +  STORE  INCREMENTS  +  DCYB*BETA 


CFY1  =  ATAB16  =  BASIC  SIDE  FORCE  COEFFICIENT  -  CY(BETA) 

EPA02  =  ATAB21  =  BETA  MULTIPLIER  TABLE 

CYDAD  =  ATAB75  =  SIDE  FORCE  COEFFICIENT  DUE  TO  AILERON  DEFLECTION 
DAILD  =  AILERON  DEFLECTION  (DEG) 

CYDRD  =  ATAB69  =  SIDE  FORCE  COEFFICIENT  DUE  TO  RUDDER  DEFLECTION 
DRUDD  =  RUDDER  DEFLECTION  (DEG) 

DRFLX5  =  ATAB88  =  FLEX  MULTIPLIER  ON  CYDRD  (=0.89) 

EPA43  =  ATAB30  =  MULTIPLIER  ON  CNDR,  CLDR,  CYDR  DUE  TO  SPEEDBRAKE 
(=1.0) 


CYDTD  =  ATAB72  =  SIDE  FORCE  COEFFICIENT  DUE  TO  DIFFERENTIAL  TAIL 
DEFLECTION  -  CYDDT 

DTFLX5  =  ATAB10  =  FLEX  MULTIPLIER  ON  CYDTD  (=0.975) 

DTFLX6  =  ATAB77  =  FLEX  INCREMENT  TO  CYDTD  (=0.0) 

DTALD  =  DIFFERENTIAL  TAIL  DEFLECTION  (DEG)  WHICH  IS 

DIRECTLY  PROPORTIONAL  TO  AILERON  DEFLECTION  AND 
IS  PRIMARILY  USED  TO  ASSIST  IN  ROLLING  THE 
F-15B  (DTALD=0 . 3*DAILD) 

CFYP  =  ATAB13  =  SIDE  FORCE  COEFFICIENT  DUE  TO  ROLL  RATE  (CYP) 

PB  =  (PEOBB*SPAN)/(2*VILWF) 

PEOBB  =  ROLL  RATE  IN  RAD/SEC  =  P 
SPAN  =  WING  SPAN  =  42.8  FEET  =  SWING 
VILWF  =  VELOCITY  IN  ET/SEC  =  VTRFPS 
CFYR  =  ATAB07  =  SIDE  FORCE  COEFFICIENT  DUE  TO  YAW  RATE  (CYR) 

RB  =  (REOBB*SPAN)/(2*VILWF) 

REOBB  =  YAW  RATE  IN  RAD/SEC  =  R 
CYRB  =  ATAB93  =  ASSYMETRIC  CY  AT  HIGH  ALPHA  (ANGLE  OF  ATTACK) 

DCYB  =0.0  THERE  IS  NO  INCREMENT  DELTA  CYB  (SIDE  FORCE) 

DUE  TO  A  2-PLACE  CANOPY  ON  THE  F15B.  THIS  IS 

BECAUSE  THE  SAME  CANOPY  IS  USED  CN  BOTH  THE 
BASELINE  F15A  AND  THE  F15B.  THE  SIDEFORCE  IS  THE 
SAME  FOR  BOTH  VERSIONS  OF  THE  F15  AND  ALREADY 
INCLUDED  IN  THE  BASIC  SIDE  FORCE  (CFY1) .  THE  TWO 
PLACE  CANOPY  IS  MOUNTED  DIFFERENTLY  HOWEVER,  SO 
THERE  IS  A  DIFFERENCE  IN  YAWING  AND  ROLLING  MCMENT. 
(SEE  DCNB  AND  DCLB) 


q* **********  cp2  ***************************************************c 
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c 

C  CFZ  =  FORCE  IN  STABILITY  MIS  Z  DIRECTION  (CL  IN  BODY  MIS) 

C  CFZ  =  CFZ1  +  CZDSPD  +  STORE  INCREMENTS  +  DCL*BETA 

C 

C 

C  CFZ1  =  ATAB17  =  BASIC  LIFT  COEFFICIENT  -  CL 
C  CZDSPD  =  ATAB26  =  DELTA  CL  DUE  TO  SPEEDBRAKE 

C  SET  TO  0  DUE  TO  THE  REASONS  GIVEN  ABOVE  IN  CXDSPD 

C  DCL  =  BTAB01  =  DELTA  CL  DUE  TO  2 -PLACE  CANOPY  (F15B)  (=0.0) 

C 

C**********  CLM  ******************************************************C 


C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


CML  =  TOTAL  ROLLING  MOMENT  COEFFICIENT  IN  BODY  MIS 
CML  =  CML1*EPA02  +  CLDAD*DAILD  +  [ CLDRD*DRUDD*DRFLX1 ] *EPA43  + 
[CLDTD*DTFLX1  +  DTFLX2]*ETALD  +  CMLP*PB  +  CMLR*RB  + 

STORE  INCREMENTS  +  CLDSPD  +  DCLB*BETA 
CML1  =  ATAB01  =  BASIC  ROLLING  MOMENT  COEFFICIENT  -  CL(BETA) 

EPA02  =  ATAB21  =  BETA  MULTIPLIER  TABLE 

CLDAD  =  ATAB73  =  ROLL  MOMENT  COEFFICIENT  DUE  TO  AILERON  DEFLECTION 
-(CLDA) 

DAILD  =  AILERON  DEFLECTION  (DEG) 

CLDRD  =  ATAB67  =  ROLLING  MOMENT  COEFFICIENT  DUE  TO  RUDDER 
DEFLECTION  -(CLD) 

DRUDD  =  RUDDER  DEFLECTION  (DEG) 

DRFLX1  =  ATAB80  =  FLEX  MULTIPLIER  ON  CLDRD  (=0.85) 

EPA43  =  ATAB30  =  MULTIPLIER  ON  CNDR,  CLDR,  CYDR  DUE  TO  SPEEDBRAKE 
(*1.0) 

CLDTD  =  ATAB70  =  ROLL  MOMENT  COEFFICIENT  DUE  TO  DIFFERENTIAL  TAIL 
DEFLECTION  -  CLDD 

DTFLX1  =  ATAB04  =  FLEX  MULTIPLIER  ON  CLDTD  (=0.975) 

DTFLX2  =  ATAB84  =  FLEX  INCREMENT  TO  CLDTD  (=0.0) 

DTALD  =  DIFFERENTIAL  TAIL  DEFLECTION  (DEG)  WHICH  IS 

DIRECTLY  PROPORTIONAL  TO  AILERON  DEFLECTION  AND 
IS  PRIMARILY  USED  TO  ASSIST  IN  ROLLING  THE  F-15B 
(DTALD  =  0.3*DAILD) 

CMLP  =  ATAB02  =  ROLL  DAMPING  DERIVATIVE  -CLP 

PB  =  (PEOBB*SPAN)/(2*VILWF) 

PEOBB  =  ROLL  RATE  IN  RAD/SEC  =  P 
SPAN  =  WING  SPAN  =42.8  FEET  =  SWING 
VILWF  =  VELOCITY  IN  FT/ SEC  =  VTRFPS 
CMLR  =  ATAB11  =  ROLLING  MOMENT  COEFFICIENT  DUE  TO  YAW  RATE  -  CLR 

RB  =  (REOBB*SPAN)/ (2*VILWF) 

REOBB  =  YAW  RATE  IN  RAD/SEC  =  R 
CLDSPD  =  ATAB29  =  DELTA  CL  DUE  TO  SPEEDBRAKE 

SET  TO  0  DUE  TO  THE  REASONS  GIVEN  ABOVE  IN  CXDSPD 
DCLB  =  BTAB04  =  INCREMENT  DELTA  CLB  (ROLLING  MOMENT)  DUE  TO  2-PLACE 
CANOPY  FROM  PSWT  499 


C***********  ***************************************************C 
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CMM  =  TOTAL  PITCHING  MCMENT  COEFFICIENT  IN  STABILITY  AXIS 
(BODY  AXIS  -  AS  WELL) 

CMM  =  CMM1  +  CMM0*QB  +  STORE  INCREMENTS  +  CMDSPD  +  DCM 


CMM1  =  ATAB03 
CMMQ  =  ATAB05 
QB 


CM 


C 
C 

c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 

C***********  cnm  ***************************************************c 

c 

CMN  =  TOTAL  YAWING  MCMENT  COEFFICIENT  IN  BODY  AXIS 
CMN  =  CMN1*EPA02  +  CNDAD*DAILD  +  [ CNDRD*DRUDD*DRFLX3 ] *EPA43 

+[CNDTD*DTLX3  +  DTFLX4]*DTALD  +  CMNP*PB  +  CMNR*RB  +  CURB 
+DCNB2*EPA36  +  STORE  INCREMENTS  +  CNDSPD  +  DCNB*BETA 
CMN1  =  ATAB12  =  BASIC  YAWING  MCMENT  COEFFICIENT  -  CN  (BETA) 


BASIC  PITCHING  MOMENT  COEFFICIENT 
PITCH  DAMPING  DERIVATIVE  -  CMQ 
(C;EOBB*MAC)  /  ( 2*VILWF) 

QIOBB  =  PITCH  RATE  IN  RAD/ SEC  =  Q 
MAC  =  MEAN  AERODYNAMIC  CHORD  =  15.94  FEET  =  OWING 
VILWF  =  VELOCITY  IN  FT/ SEC  =  VTRFPS 
CMDSPD  =  ATAB25  =  DELTA  CM  DUE  TO  SPEEDBRAKE 

SET  TO  0  DUE  THE  REASONS  GIVEN  ABOVE  IN  CXDSPD 
DCM  =  BTAB02  =  DELTA  CM  DUE  TO  2-PLACE  CANOPY  (F15B)  (=0.0) 


DAILD 

CNDRD 

DRUDD 


EPA02  =  ATAB21  =  BETA  MULTIPLIER  TABLE 
CNDAD  =  ATAB74  =  YAW  MOMENT  COEFFICIENT  DUE  TO  AILERON 
DEFLECTION  -CNDA 

=  =  AILERON  DEFLECTION  (DEG) 

=  ATAB68  =  YAWING  MCMENT  COEFFICIENT  DUE  TO  RUDDER 
DEFLECTION  -CNDR 
=  RUDDER  DEFLECTION  (DEG) 

DRFLX3  =  ATAB85  =  FLEX  MULTIPLIER  ON  CNDRD 

EPA43  =  ATAB30  =  MULTIPLIER  OS  CNDR,  CLDR,  CYDR  DUE  TO  SPEEDBRAKE 
CNDTD  =  ATAB71  =  YAWING  MCMENT  COEFFICIENT  DUE  TO  DIFFERENTIAL  TAIL 
DEFLECTION  -  CNDDT 

DTFLX3  =  ATAB08  =  FLEX  MULTIPLIER  ON  CNDTD 
DTFLX4  =  ATAB09  =  FLEX  INCREMENT  CN  CNDTD  (=0.0) 

=  DIFFERENTIAL  TAIL  DEFLECTION  (DEG)  WHICH  IS 
DIRECTLY  PROPORTIONAL  TO  AILERON  DEFLECTION 
AND  IS  PRIMARILY  USED  TO  ASSIST  IN  ROLLING 
THE  F-15B  (DTALD  =  0.3*DAILD) 

YAWING  MCMENT  COEFFICIENT  DUE  TO  ROLL  RATE  -  CNP 
( PEOBB*SPAN ) / ( 2*V ILWF ) 

PEOBB=ROLL  RATE  IN  RAD/ SEC  =  P 
SPAN  =  WING  SPAN  =  42.8  FT  =  SWING 
VILWF  =  VELOCITY  IN  FT/ SEC  =  VTRFPS 
YAW  DAMPING  DERIVATIVE  -  CNR 
(RE0BB*SPAN)/(2*VILWF) 

REOBB  =  YAW  RATE  IN  RAD/SEC  =  R 
CNRB  =  ATAB86  =  ASSYMETRIC  CN  AT  HIGH  ALPHA 

DCNB2  =  ATAB44  =  DELTA  CNB  WITH  STABILATOR  EFFECT  -  DELCNB  (=0.0) 
EPA36  =  ATAB94  =  MULTIPLIER  ON  DCNB2  (=BETA) 

CNDSPD  =  ATAB28  =  DELTA  CN  CUE  TO  SPEEDBRAKE 

SET  TO  0  DUE  TO  THE  REASONS  GIVEN  ABOVE  IN  CXDSPD 
DCNB  =  BTAB05  =  INCREMENT  DELTA  CNB  (YAWING  MCMENT)  DUE  TO 


DTALD  = 


CMNP  =  ATAB06 
PB 


CMNR  =  ATAB14 
RB 
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2-PLACE  CANOPY  (F15B) 

C*******************************************************3**************C 

C 

C  MISCELLANEOUS  COEFFICIENTS  AND  NAME  CONVERSION  TABLE 
C 


c 

1988  F15 

ORIGINAL 

c 

AEROBASE  NAME 

PROGRAM  NAME 

DEFINITION 

c 

************* 

************ 

********** 

c 

c 

AL77D 

AL 

ANGLE  OF  ATTACK 

c 

(DEG) 

c 

BE77D 

BETA 

SIDESLIP  ANGLE 

c 

(DEG) 

c 

BE77D 

RBETA 

SIDESLIP  ANGLE 

c 

(RAD) 

c 

B077D 

ABET 

ABSOLUTE  VALUE  OF 

c 

SIDESLIP  ANGLE 

c 

(DEG) 

c 

DAILA 

DAILA 

ABSOLUTE  VALUE  OF 

c 

AILERON  DEFLEC¬ 

c 

TION  (DEG) 

c 

DAILD 

DDA 

AILERON  DEFLEC¬ 

c 

TION  (DEG) 

c 

DRUABS 

ARUD 

ABSOLUTE  VALUE  OF 

c 

RUDDER  DEFLEC¬ 

c 

TION  (DEG) 

c 

DRUABS 

RARUD 

ABSOLUTE  VALUE  OF 

c 

RUDDER  DEFLEC¬ 

c 

TION  (RAD) 

c 

DRUDD 

DRUDD 

RUDDER  DEFLECTION 

c 

(DEG) 

c 

DSTBD 

DELESD(R) 

AVERAGE 

c 

STABILATOR 

c 

DEFLECTION 

c 

DEG  (RAD) 

c 

DTALD 

DELECD(R) 

DIFFERENTIAL  TAIL 

c 

DEFLECTION 

c 

DEG  (RAD) 

RBETA=BETA/ DEGRAD 
DAILA=ABS(DDA) 

C 


PB=(P*BWING ) / ( 2 . OdO*VTRFPS) 

QB=(Q*CWING)/(2.0d0*VTRFPS) 

RB= (R*BWING) / ( 2 . OdO*VTRFPS) 

C 

C  THE  F-15B  AERO  DATA  TABLES  DO  NOT  CONTAIN  STABILITY  COEFFICIENT 
C  DATA  FOR  BETA  AND  RUDDER  DEFLECTION  ,DRUDD,  LESS  THAN  0 

C  DEGREES.  THE  ABSOLUTE  VALUE  OF  BETA,  ABET,  AND  THE  ABSOLUTE 

C  VALUE  OF  RUDDER  DEFLECTION,  ARUDD,  ARE  USED  IN  THE  FOLLOWING 
C  EQUATIONS.  IN  RADIANS  THESE  PARAMETERS  ARE  RABET  AND  RARUD, 

C  RESPECTIVELY.  IN  SOME  CASES  THE  COEFFICIENT  IS  MULTIPLIED  BY  A 

C  -1  FOR  PARAMETER  VALUES  LESS  THAN  ZERO. 
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c 

C  EPA02  IS  A  MULTIPLIER  THAT  ADJUSTS  THE  PARTICULAR  COEFFICIENT 
C  IT  IS  WORKING  ON  (CFY1,CML1,CMN1)  BY  CHANGING  THAT  PARTICULAR 
C  COEFFICIENTS  SIGN  (POSITIVE  OR  NEGATIVE)  DEPENDENT  ON  THE  SIGN 
C  OF  THE  SIDESLIP  ANGLE  (BETA) .  IF  BETA  IS  NEGATIVE  THEN 
C  EPA02=-1.0.  IF  BETA  IS  POSITIVE  THEN  EPA02=1.0.  SINCE  THIS 
C  FUNCTION  IS  DISCONTINUOUS  AT  TOE  ORIGIN  A  CUBIC  SPLINE  HAS 
C  BEEN  EMPLOYED  TO  REPRESENT  THIS  EUNCTION  IN  ORDER  THAT 
C  AUTO  CAN  RUN. 

C 

IF  (BETA  .LT.  -1.0)  THEN 
EPA02S=  -l.OdO 
ENDIF 
C 

IF  ((BETA  .GE.  -1.0)  .AND.  (BETA  .LE.  1.0))  THEN 
EPA02S=-1 . 0d0+ ( 1 . 50d0* ( (BEIA+1 . OdO ) **2 . OdO ) ) - 
1  (0.50d0*UBETA+1.0d0)**3.0d0)) 

ENDIF 

C 

IF  (BETA  .GT.  1.0)  THEN 

EPA02S=1.0d0 

ENDIF 

C 

IF  (BETA  .LT.  -5.0)  THEN 
EPA02L=  -l.OdO 
ENDIF 
C 

IF  ((BETA  .GE.  -5.0)  .AND.  (BETA  .LE.  5.0))  THEN 
EPA02L=-1 . 0d0+ (0 . 060d0* ( (BETA+5 . OdO ) **2 . OdO ) ) - 
1  ( 0 . 0040d0* ( ( BETA+5 . OdO ) **3 . OdO ) ) 

ENDIF 

C 

IF  (BETA  .GT.  5.0)  THEN 
EPA02L=1 . OdO 
ENDIF 

*****  DIFFERENTIAL  ELEVATOR  ***** 

DTALD=0 . 30d0*DAILD 
DELEDD=0 . 30d0*DDA 
DELEDR=0 . 30d0*  (DDA/DEX3RAD) 

C 

<2**********  cp2  *************************************************** c 
C 

CFZl=-0. 0036937 6+( 3. 78028702*RAL)+(0.6921459*RAL*RAL)-( 5. 0005867 
+*(RAL**3))+(1. 94478199* (RAL**4) )+(0 . 40781955*DELESR)+(0 . 1U114579 
+*(DELESR*DELESR) ) 

C 

CFZ=CFZ1 

C 

C*********  CFX  ******«*********************************************c 
C 

CL=CFZl/57. 29578 
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non  ooooooo 


THIS  CONVERSION  OF  CFZ1  TO  CL  IS  AN  ARTIFACT  FROM  THE 
CURVE  FITTING  PROCESS  WHERE  ALL  THE  INDEPENDENT  VARIABLES 
WERE  ANGLES  THAT  WERE  CONVERTED  FROM  DEGREES  TO  RADIANS. 

IT  JUST  SO  HAPPENED  THAT  FOR  CFX1  ONE  OF  THE  VARIABLES 
WAS  NOT  AN  ANGLE  BUT  A  DIMENSIONLESS  COEFFICIENT. 

cm=0 .0180682! -K0 .01556573*CL)+(  498 . 96208868*CL*CL) 

+-(14451. 56518396* (CL**3) )+{ 2132344. 6184755* (CL**4) ) 

TRANSITIONING  FRCM  LOW  AQA  DRAG  TABLE  TO  HIGH  ADA  DRAG  TABLE 

CFX2=0. 0267297- (0.1Q646919*RAL)+ (5. 39836337*RAL*RAL) 

+-( 5. 0086893* (RAL**3) )+(1.34148193*(RAL**4) )+ 

+(0 . 20978902*DELESR)+(0 . 3CS04211*(DELESR**2) )+0 . 09833617 
C 

Al=20 . OdO/DEGRAD 
A2=30 . OdO/DEGRAD 
A12=A1+A2 

BA=2 . 0/ ( -A1**3t3 . *A1*A2* (A1-A2 ) +A2**3 ) 

BB=~3. OdO*BA* ( A1+A2 ) / 2 . OdO 
BC=3 . 0d0*BA*Al*A2 
BD=BA*A2**2* ( A2-3 . OdO*Al ) /2 . OdO 
F1=BA*RAL**3+BB*RAL**2+BC*RAL+BD 
F2=-BA*RAL**3+ ( 3 . 0d0*A12*BA+BB)*RAL**2 . OdO- 
+  (BC+2 . 0d0*A12*BB+3 . 0d0*A12**2*BA)*RAL+ 

+  BD+A12*BC+A12**2*BB+A12**3*BA 

C 

IF  (RAL  .LT.  Al)  THEN 

C 

CFX=CFX1 

C 

ELSEIF  (RAL  .GT.  A2)  THEN 
C 

CFX=CFX2 

C 

ELSE 

C 

CFX=CFX1*F1+CFX2*F2 

C 

ENDIF 

C 

C*********  cpy  *********************x*****************************c 

C 

DTFLX5-0 . 975d0 
DRFLX5=0 . 89d0 
C 

CFYl=-0. 050 60386- (0.12342073*RAL)+(1.04501136*RAL*RAL) 

+-(0 .17239516* (RAL**3) )-( 2 . 90979277*(RAL**4) ) 
++(3.06782935*(RAL**5))-(0.8S422116*(RAL**6)) 

+- (0 . 06578812*RAL*RABET) -( 0 . 71521988*RABET) -( 0 . 00000475273 
+*(RABET**2) ) - ( 0 . 04856168*RAL*DELESR) - ( 0 . 05943607*RABET*DELESR ) + 
+ ( 0 . 02018534*DELESR) 
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c 

IF  (RAL  .LT.  .52359598)  THEM 
C 

CFYP=0 . 01460 6188+ (2 . 52405055*RAL) - ( 5 . 02687473* (RAL**2 ) ) 
+-(106.43222962*(RAL**3))+(256.80215423*(RAL**4)) 

++(1256. 39636248* (RAL**5)) 

+•■  ( 3887 . 92878173*(RAL**6) ) - f 2863 . 16083460* (RAL**7 ) ) + 

+(17382. 72226362*(RAL**8))-(13731.65408408*(RAL**9)) 

ENDIF 

C 

IF  ((RAL  .GE.  .52359998)  .AND.  (RAL  .LE,  .610865))  THEN 

C 

CEYP=0 . 00236511+ ( 0 . 52044678* (RAL-0 . 52359998) )-(12 . 8597002* (RAL- 
+0 . 52359998 ) **2 ) +( 75 . 46138* (RAL-0 . 52359998 )**3 ) 

ENDIS 

C 

IF  (RAL  .GT.  0.610865)  THEN 
C 

CFYP=0.0d0 

ENDIF 

C 

IF  (RAL  .LT.  -0.06981)  THEN 
C 

CF7P=0 . 35d0 
ENDIF 
C 

IF  ((RAL  .GE.  -0.06981)  .AND.  (RAL  .LT.  0.0))  THEN 
C 

CF?R=0 . 34999999+ ( 35 . 4012413* (RALrO . 0 6981 ) **2 ) - ( 493 . 33441162* 
+(RAL+0. 06981 )**3) 

ENDIF 

C 

IF  ((RAL  .GE.  0.0)  .AND.  (RAL  .LE.  0.523599))  THEN 

C 

CFYR=0 . 35468605- ( 2 . 26998141*RAL) + ( 51 . 82178387*RAT*RAL ) 
+-(718.55069823*(RAL**3)) 

++( 4570 . 00492172*  (RAL**4 ))-( 14471 . 38028351*  (RAL*  ■'•5 )) + 

+(22026 . 58930662* (RaL**6) ) - ( 12795 . 99029404*(RAL**7 ) ) 

ENDIF 

C 

IF  ((RAL  .GT.  0.523599)  .AND.  (RAL  ,il2.  0.61087))  THEN 
C 

CF7R=0 . 00193787+(1 .78332495* (RAL-0 . 52359903) )-(41 . 63198853* (RAL- 
+0 . 52359903)**2 )+( 239 . 97909546* (RAL-0 . 52359903)**?) 

ENDIF 

C 

IF  (RAL  .GT.  0.61087)  THEN 

CE7R=0.0d0 

ENDIF 

C 

IF  (RAL  .LT.  0.55851)  THEN 

n 


88 


c 

c 


c 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 


CYDAD=-0.00020812+(0.0Q062122*RAL)+(0.00260729*RAL*RAL) 

++ (0 . 00745739* (RAL**3 ) ) - ( 0 . 0365611* (RAL**4 ) ) 

+- ( 0 . 04532683* (RAL**5 ) ) + (0 . 2067 4845* (RAL**6) ) 

+- ( 0 . 13264434* (RAL**7 ) ) - ( 0 . 00193383* (RAL**8) ) 

ENDIF 

IF  ((RAL  .GE.  0.55851)  .AND.  (RAL  .LT.  0.61087))  THEN 

CYDAD=0.00023894+(0. 00195121*(RAL-0. 55851001) )+(0. 02459273 
+* (RAL-0 . 55851001 ) **2 ) - ( 0 . 1202244* ( (RAL-0 . 55851001)**3) ) 

ENDIF 

IF  (RAL  .GE.  0.61087)  THEN 

CYDAD=0 . 27681285- (2 . 02305395*RAL)+( 6 . 01180715*RAL*RAL) 

+- ( 9 . 24292188* (RAL**3) )+(7 . 59857819*(RAL**4) ) 

+-(2 . 8565527* (RAL**5) )+(0 . 25460503* (RAL**7 ) ) 

+- ( 0 . 01819815* (RAL**9) ) 

ENDIF 


IF  (RAL  .LE.  0.0)  THEN 

EPA43=1.0d0 

ENDIF 

IF  (RAL  .GT.  0.0  AND  .LE.  0.6283185)  THEN 
0.6283185  RADIANS  =  36  DEGREES 

EPA43=0 . 9584809+ (4 . 13369452*RAL) - ( 18 . 31288396*RAL*RAL) + 

+ ( 19 . 5511466* (RAL**3 ) ) - ( 1 . 09295946*RAL*DSPBD) + ( 0 . 17441033* 
+DSPBD*DSPBD ) 

ENDIF 

IF  (RAL  .GT.  0.6283185)  THEN 

EPA43=1.0d0 

ENDIF 

********************************************************* 

*  NOTE  -  THE  PARAMETER  EPA43  IS  A  MULTIPLIER  ON  RUDDER  * 

*  EFFECTIVENESS  DUE  TO  SPEEDBRAKE.  THIS  TABLE  IS  ALSO  * 

*  LIMITED  TO  36  DEG  AO  A.  HOWEVER,  THERE  IS  NO  AERODY  * 

*  NAMIC  EFFECT  FOR  ANGLES  OF  ATTACK  LESS  THAN  16  DEG,  * 

*  AND  THE  SPEEDBRAKE  IS  AUTOMATICALLY  RETRACTED  AT  ADA  * 

*  GREATER  THAN  15  DEG.  THEREFORE,  THIS  TABLE  SHOULD  * 

*  NOT  BE  NECESSARY  FOR  THE  ORDINARY  OPERATION  OF  THE  * 

*  AIRCRAFT  * 

**********  *********************************************** 


CYDRD=0.00310199+(0.00119963*RAL)+(0.02806933*RAL*RAL) 

+-(0.12408447*(RAL**3))-(0.12032121*(RAL**4)) 

++(0. 7915027 9*(RAL**5))-(0.86544347*(RAL**6)) 

++(0 . 27845115* (RAL**7 ) )+(0 .00122999*RAL*RARUD)+(0 .00145943 
+*RARUD) - ( 0 . 01211427*RARUD*RARUD)+ ( 0 . 00977 937* ( RARUD**3 ) ) 

CYETD=-0.00157745-(0.0020881*RAL)+(0.00557239*RAL*RAL) 

+- ( 0 . 00139886* (RAL**3 ) )+( 0 . 04956247* (RAL**4 ) ) 

+- (0 . 0135353*(RAL**5))-(0.11552397*(RAL**6) ) 
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++ ( 0 . 11443452* (RAL**7 ) ) - (0 . 03072189* (RAL**8 ) ) - (0 . 01061113* 

+ (RAL**3 ) *DELESR)- ( 0 . 00010529*RAL*RAL*DELESR*DELESR) 

«~(0 . 00572463*RAL*DELESR*DELESR) 

++  ( 0 . 01885361*RAL*RAL*DELESR) -  ( 0 . 01412258*RAL* (DELESR** ?  ) ) 

+-  (0 . 00082.7  76*DELESR) +(0 . 00404354*  (DELESR**2 ) )  - 
+(0 .00212189*(DELESR**3) )+(G.00655063*fDELESR**4) ) 

++( 0 ..  03341584*  (DELESR**5  ) ) 

C 

RALY1=0 . 6108652 
RALY2=90 . OdO/DEGRRD 
FBETYl=-0. 0872565 
RBETY2=0. 17 45329 
C 

AY=0 . 1640d0 
A£?TARY=0. 95993 
BSTARY=0. 087266 
C 

ZEmy=  ( 2 . 0D0*ASTAKY- (RALY1+RALY2 ) )/  (RALY2-RALY1 ) 

ETAY=  ( 2 .  ODO*BSTr;  SCI-  ( RBETY1+RBETY2 ))/( RBETY2-RBETY1 ) 

C 

X= ( 2 . ODO*RAL- (RALY1+RALY2 ) )/(RALY2-RALYl ) 

Y=(2. ODO*RBETA- (RBETY1+RBETY2 ) ) / ( RBETY2-RBETY1 ) 

C 

FY=(  (5.  ODO*  ( ZETAY**2  )  )  -  (  4 . 0D0*ZE3AY*X)  -1 .  ODO  )*  ( (  (X**2 -I .  ODO ) 
+**2 ) * ( 1 .  ODO/  ( ( ( ZETAY**2 )  -1 .  ODO  )  **3 ) 

C 

GY=  (( 5 .  ODO*  (ETAY**2  ))-(  4 .  ODO*ETKY*Y) -1 .  ODO )-'(( (Y**2 ) -1 .  CDO  )  **2 ) 
+* ( 1 . ODO/ ( ( ( EIAY**2 ) -1 . ODO ) **3 ) ) 

-t 

CYRB=AY*F7*G¥ 

C 

IF  (RAL  .LT.  0.6108652)  TFjiN 
C 

CYRB=0.0d0 
GOTO  500 
ENDIF 
C 

IF  ( (RBETA  .LT.  -0.0872665)  .OR.  (RBETA  .GT.  0.1745329))  THEN 
C 

CYR;  O.OdO 
GOTO  500 
ENDIF 
C 

jOO  CFY=  ( CFY1*EIPA0  2l)+(  CYDAD*Dr  A ) + ( CYDRD*DRUDD*DRFLX5*EPA43 ) + 

+( (Cn7TD*DTFLX5)'lfDELE23D)+(CEYP*PB)+(CFra*RB) 
t+CYRB 
C 

C************  cm  ***************^-********************************c 

c 

DTFLXl=0.9750d0 

DFFLXl=0.850d0 

C 

CMLl=-0 . 00238235- ( 0 . 046 \6235*RAL)+ ( 0 . 1Q553168*RAL*RAL ) 
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++ ( 0 . 10541585* (RAL**3 ) ) - ( 0 . 402547 65* (RAL**4 ) ) 

++(0 . 32530491*(RAL**5 ) )-(0 . 08496121*(RAL**6) ) 

++(0.00112288*(RAL**7))-(0.05940477*RABE3*RAL)- 

+(0 .07356236*RABETT)-(0 .00550119*RABET*RABEn,)+(0 . 00326191 

+*(RABET**3)) 

C 

IF  (HAL  .LT.  0.29671)  THEN 
C 

CMLP=-0 . 24963201-(0. 03106297*RAL)+(0 . 12430631*RAL*RAL) 

+- (8. 95274618* (RAL**3) )+(100 . 33109929*(RAL**4) ) 

++(275 . 70069578* (RAL**5) )-(1178 . 83425699*(RAL**6) ) 

+-(2102 . 66811522*(RAL**7) )+( 2274. 89785551* (RAL**8) ) 

ENDIF 

C 

IF  ((RAL  .GE.  0.29671)  .AND.  (RAL  .LT.  0.34907))  THEN 

C 

CMLP=-0 . 1635261- (3 . 77847099* (RAL-0 . 29671001) )+(147 . 47639465 
+* (RAL-0 . 29671001 ) **2 ) - ( 1295 . 947 99805* (RAL-0 . 29671001 ) **3 ) 

ENDIF 

C 

IF  (RAL  .GE.  0.34907)  THEN 

C 

CMLP=-1 . 37120291+(7 . 0611218.''  *RAL)-  (13 . 57010422*RAL*RAL) 

++( 11 . 21323850* (RAL**3 ) ) 

+- ( 4 . 26789425*(RAL**4) )+(0 . 6237381* (RAL**5) ) 

ENDIF 

C 

IF  (RAL  .LT.  0.7854)  THEN 
C 

CMLR=0 . 03515391+ ( 0 . 59296381*RAL) + ( 2 . 27456302*RAL*RAL ) 

+- (3. 8097803* (RAL**3)) 

+- ( 45 . 83162842* (RAL**4 ) ) + ( 55 . 31669213* (RAL**5 ) ) + 

+(194. 29237485*(RAL**6))-( 393. 22969953*(RAL**7))+(192. 20860739* 
+(RAL**8) ) 

ENDIF 

C 

IF  ((RAL  .GE.  0.7854)  .AND.  (RAL  .LE.  0.87266))  THEN 
C 

CMLR=0 . 0 92557 9071- ( 0 . 6000000238* (RAL-0 . 7853999734) ) 
++(1.3515939713*( (RAL-0. 7853999734)**2)) 

++(29. 0733299255*( (RAL-0. 7853999734)**3)) 

ENDIF 

C 

IF  (RAL  .GT.  0.87266)  THEN 
C 

CMLR= -311 . 126041+ ( 1457 . 23391042*RAL) - ( 2680 . 19461944*RAL*RAL)  + 
K  2361 . 44914738* ( RAL**3 )) - ( 893 . 63567263* ( RAL**4 ))+( 68 . 23501924* 
+(RAL**6))-(1.72572994*(RAL**9)) 

ENDIF 

C 
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CLDAD=0 . 00057 626+ (0 . 00038 47 9*RAL ) - ( 0 . 00502091*RAL*RAL) 
++(0.00161407*(RAL**3))+(0.02268829*(RAL**4)) 

+- (0. 03935269* (RAL**5))+(0.02472827*(RAL**6) ) 

+- (0 . 00543345*(RAL**7 ) ) +(0 . 0000007520348*DELESR*RAL ) + 

+(0 .000000390773*DELESR) 

C 

CLDRD=0 . 00013713- (0 . 00035439*RAL) -(0 . 00227912*RAL*RAL) 

++(0 . 00742636* (RAL**3) )+(0 .00 991839* (RAL**4) ) 
+-(0.04711846*(RAL**5))+(0.046124*(RAL**6)) 

+- ( 0 . 0137 9021* (RAL**7 ) ) + (0 . 00003678685*RAKUD*RAL ) + 
+(0.00001043751*RARUD)-(0.00015866*RAKUD*RARUD)+(0. 00016133 
+*(RARUD**3)) 

C 

CLDTD=0 . 00066663+{0. 00074174*RAL)+(0 . 00285735*RAL*RAL) 

+- ( 0 . 02030692* (RAL**3 ) ) - ( 0 . 00352997* (RAL**4) ) 
++(0.0997962*(RAL**5))-(0. 14591227* 

+(RAL**6) )+(0.08282004*(RAL**7) ) 

+-(0.0168667*(RAL**8) )+(0. 003061 42*(RAL**3)*DELESR) 

+- ( 0 . 00110266*RAL*RAL* ( DELESR**2) ) + ( 0 . 00088031*RAL* 

+(DELESR**2 ) ) - (0 . 00432594*RAL*RAL*DELESR) - 
+(0 .00720141*RAL*(DELESR**3) ) 

+- ( 0 . 00034325*DELESR) + (0 . 00033433* (DELESR**2 ) )+(0 . 00800183 
+* ( DELESR**3 )) - ( 0 . 00555986* ( DELESR**4 ))-(0. 0184117 2*( DELESR**  5 ) ) 

C 

IF  (RAL  .LT.  0.0)  THEN 

C 

DCLB=-0 . 000060d0 

ENDIF 

C 

IF  ((RAL  .GE.  0.0)  .AND.  (RAL  .LE.  0.209434))  THEN 
C 

DCLB=-0 . 000060d0+(0 . 0041035078*RAL*RAL)-(0 . 0130618699*(RAL**3) ) 

ENDIF 

C 

IF  (RAL  .(7T.  0.209434)  THEN 
C 

DCLB=0 . CdO 

ENDIF 

C 

CML=  ( CML1*EPA0  2  S ) + (  CLDAD*DDA) + ( CLERD*DRUDD*DRFLX1*EPA43 ) + 

+  ( ( CLDTD*DTFLX1 ) *DELEDD ) + ( CMLP' *  PB ) + ( CMU**RB ) + ( DCLB*BETA ) 

C 

(^**************  *★*******★*********★**************************£ 

C 

CMM1=0 . 00501496-(0. 08004901*RAL)- (1 . 03486675*RAL*RAL) 

+-(0 . 68580677* (RAL**3) )+( 6. 46858488*(RAL**4) ) 

+-(10 . 15574108* (RAL**5) )+ 

+( 6. 44350808*(RAL**6) )-(l. 46175188*(RAL**7) ) 

++(0 . 24050902*RAL*DELESR) 

+- ( 0 . 42629958*DELESR) - (0 . 03337 449*DELESR*DELESR) 

+-(0 . 53951733*(DELESR**3) ) 

modified  25  Jul  90  to  use  new  curve  fit  for  CMQ 
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c 

C  OLD  EQUATION 

C 

C  IF  (RAL  .LE.  0.25307)  THEN 

C 

C  CMMQ=-3 . 8386262+(13 . 54661297*RAL)+(402 . 5301155 9*RAL*RAL ) 

C  +- ( 6660 . 95327122* (RAL**3 ) ) - ( 62257 . 89908743* (RAL**4 ) ) 

C  ++(261526. 10242329*(RAL**5)) 

C  ++(2177190 . 33155227* (RAL**6 ) ) - (703575 . 13709062*(RAL**7 ) ) - 
C  +(20725000 . 34643054* (RAL**8) )- ( 27829700 . 5333364?*(RAL**9) ) 

C  ENDIF 

C 

C  IF  ((RAL  .GT.  0.25307)  .AND.  (RAL  .LT.  0.29671))  THEN 
C 

C  CMMQ=-8 . 4926528931- ( 2705 . 3000488281*(RAL-0 . 2530699968 ) ) 

C  ++(123801 . 5*(RAL-0 . 2530699968)**2) 

C  +- ( 1414377* (RAL-0 . 2530699968 ) **3 ) 

C  ENDIF 

C 

C  IF  (RAL  .GE.  .29671)  THEN 

C 

C  CMMQ=47.24676075-(709. 60757056*RAL)+(3359.08807193*RAL*RAL)- 

C  + ( 7565 . 32017266*(RAL**3 ) ) +( 8695 . 1858091* (RAL**4 ) ) 

C  +- ( 4891 . 77183313* (RAL**5 ) ) + ( 1061 . 55915089* (RAL**6 ) ) 

C  ENDIF 

C 

C  CMMQ  vs.  alpha  n  degrees 

C 

C  NEW  EQUATION 

C 

C  convert  alpha  to  degrees 

C 

A=RAL*DEGF7\D 

C 

Fl=-4 . 33509d0+A*( -0 . 141624dO+A* ( 0 . 0946448d0+A*( -0 . 00798481d0 
+  +A*(-0 ,00168344dO+A*(0 .000260037d0+A*( 6. 64054d-6+A*( 

+  -2.20055d-6+A*(-2.74413d-8+A*(7.14476d-9+A* 

+  2 . 07046d-10 ))))))))) 
c 

F2=-302 . 567+a*( 106 . 288+a*( -14 . 7034+A* ( 1 . 02524+A*( -0 . 0393491 
+  +A*  ( 0 . 00084082+A* ( -9 . 365e-6+A*4 . 2355e-8 ) ) ) ) ) ) 

F3=1724 . 99+A*( -158 . 944+A*(5 . 59729+A*(-0 .0949624+A*( 

+  0 . 000779066+A*(-2 . 47982e-6) ) ) ) ) 
c 

c  ramp  functions 
c 

Rl=l . 0-0 . 75* ( A-10 . 0 ) **2+0 . 25* ( A-10 . 0 ) **3 
R2=1.0-R1 

R3=l . 0-7 . 5* (A- 40 . 0)**2/62 . 5+(A-40 . 0 )**3/62 . 5 
R4=l . 0-R3 
c 

IF(A.LT.10.0),IHEN 
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CMMQ=F1 

ELSEIF(  A .  LT .  12 . 0  )7HEN 
OHQ=Fl*Rl+F2*R2 
ELSEIF(A. LT. 40 . 0)1HEH 
CMMQ=F2 

ELSEIF(  A.  LT .  45 . 0  )THEN 
CMMQ=F2*R3+F3*R4 
ELSE 
CJfiQ=F3 
EUDIF 
C 

Q®4=Ofll+(CM4Q*QB) 

C 

C************  ************************************************(2 

C 

DTFLX3=0 . 9750d0 
DEFLX3=0.890dO 
C 

CMN1=0 . 01441512+ ( 0 . 02242944*FAL )  -  (  0 . 30472558* (RAL**2 ) ) 

++(0 . 14475549* (RAL**3) ) 

++ ( 0 . 93140112* (RAL**4 ) ) - ( 1 . 52168677*  (RAL**5 ) ) + 

+(0 . 90743413* (RAL**6) )-(0.16510989*(RAL**7) ) 

+- ( 0 . 0461968* (RAL**8 ) ) 

++(0 .01754292*(RAL**9) )~(0 .17553807*RAL*RABET)+ 

+ ( 0 . 15415649*RAL*RABET*DELESR) 

++(0 . 14829547* (RAL**2 ) * (RABET**2 ) ) 

+~ (0 . 11605031*(RAL**2)*RABET*DELESR) 

+-(0.06290678*(RAL**2)*(BELESR**2)) 

+-(0.01404857*(RAL**2)*(DELESR**2)) 

++ ( 0 . 07225609*RABET)-(0 . 08567087* (RABET**2 ) ) 

++ ( 0 . 0118467 4* (RABET**3 ) ) 

+-  ( 0 . 00519152*RAL*BELESR) + (0 . 03865177*RABET*raLESR) 

+ f(0 . 00062918*DELESR) 

C 

CNDRD=-0 . 00153402+(0 . 00184982*RAL)-(0 . 0068693*RAL*RAL) 

++(0 . 01772037*(RAL**3) ) 

++ ( 0 . 03263787* (RAL**4 ) ) - ( 0 . 15157163* (RAL**5 ) ) + (0 . 18562888 
+* (RAL**6) ) - ( 0 . 09661 63*(RAL**7 ) ) + ( 0 . 01859168* (RAL**8 ) ) + ( 0 . 0002587 
+*RAL*DELESR)  -  ( 0 . 00018546*RAL*DELESR*RBETA)  -  (0 . 00000517304*RBETA) 
+-  ( 0 . 00102718*RAL*RBETA)  -  ( 0 . 000068937 9*RBm*DELESR)  -(0.00040536 
+*RBETA*RAW3D)  -  ( 0 . 00000480484*DELESR*RARUD) 

+- ( 0 . 00041786*RAL*RARDD) 

++(0.0000461872*R8EEA)+(0.00434094*(RBETA**2)) 

+- ( 0 . 00490777*(RBETA**3 ) ) 

++ ( 0 . 000005157867*RARUD)+ ( 0 . 00225169*RAR0D*RARDD) - ( 0 . 00208072 
+*(RAKUD**3)) 

C 

IF  (RAL  .LT.  0.55851)  THEN 

C 

CMNP=-0 . 00635409- (1 . 14153932*RAL)+( 2 . 82119027*(RAL**2 ) )+ 

+ ( 54 . 473957 9* (RAL**3 ) ) - ( 140 . 89527 667* (RAL**4 ) ) - ( 676 . 737 46128* 
+(RAL**5) )+( 2059. 18263976*(RAL**6))+(1579. 416647 48*(RAL**7)) 
+-(8933 . 08535712*(RAL**8 ) ) +( 6806 . 54761267* (RAL**9) ) 
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ENDIF 

C 

IF  ((RAL  .GE.  0.55851001)  .AND.  (RAL  .LE.  0.61087))  THEN 
C 

CMNP=-.07023239+(1.085815*(RAL  -0.55851)) 

++( 8 . 852651* ( (RAL- . 55851) **2 ) )- (192 . 6093*( (RAL-0 . 55851)**3) ) 
ENDIF 
C 

IF  (RAL  .GT.  0.61087)  THEN 
C 

CMNP=-71 . 03693533+ (491 . 3250  6715*RAL) 

+-( 1388 . 11177  97  9* (RAL**2 ) ) + 

+(2033 . 48621905* (RAL**3 ) ) 

+-(1590 . 91322362* (RAL**4) )+(567 . 38432316*(RAL**5) ) 

+-( 44. 97702536*(RAL**7 )) +(2.8140669* (RAL**9) ) 

ENDIF 

C 

IF  (RAL  .LE.  -.069813)  THEN 
C 

CMNR=  -0 . 28050d0 
ENDIF 
C 

IF  ((RAL  .GT.  -.069813)  .AND.  (RAL  .LT.  0.0))  THEN 

C 

CMNR=-0 . 2804999948+( 35 . 9903717041* (RAL+ . 0698129982 )**2 ) 

+- ( 516 . 1574707031* (RAL+ . 0698129982 ) **3 ) 

ENDIF 

C 

IF  ((RAL  .GE.  0.0)  .AND.  (RAL  .LE.  0.78539801))  THEN 

C 

CMNR=- . 28071511-(2. 52183924*RAL)+(68 . 90860031* (RAL**2) ) 

+- ( 573 . 23100511*(RAL**3 ) )+( 2009 . 08725005* (RAL**4) ) 

+- ( 3385 . 15675307* (RAL**5 ) ) 

++(2730 . 49473149* (RAL**6) )-(848 . 12322034* (RAL**7 ) ) 

ENDIF 

C 

IF  ((RAL  .GT.  0.78539801)  .AND.  (RAL  .LT.  0.95993102))  THEN 
C 

CMNR=-0 . 1096954+ (0 . 52893072* (RAL-0 . 78539801 ) ) - ( 6 . 09109497* (RAL- 
+0 . 78539801 ) **2 ) + ( 17 . 47834015* (RAL-0 . 78539801 ) **3 ) 

ENDIF 

C 

IF  (RAL  .GE.  0.95993102)  THEN 
C 

CMNR=-0 .  HOdO 
ENDIF 
C 

CNDTD=0.00058286+(0.0007341*RAL)-(0.00746113*RAL*RAL) 
+-(0.00685223*(RAL**3) ) 

++ (0.0327 7 271* (RAL**4) )-(0.02791456*(RAL**5) ) 
++(0.00732915*(RAL**6)) 

++(0.00120456*RAL*DELESR)-(0.00168102*DELESR)+(0. 0006462* 
+DELESR*DELESR ) 
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c 

CNDAD=0 . 00008228887- (0 . 00014015*RAL) - (0 . 0013493*RAL*RAL)+ 

+(0 .00020487*(RAL**3) )+(0 .005612 41* (RAL**4) ) 
+-(0.00634392*(RAL**5)) 

++(0 . 00193323* (RAL**6) )-(2 .05815E-17*(RAL*DAILA) )+(3 .794816E-17* 
+(DAILA**3)) 

C 

DCNB=-2.500E-4 

C 

RALN1=0. 69813 
RALN2= 90 . OdO/DEGRftD 
RBETNl=-0. 174532 
RBETN2=0. 34906 
C 

AN=0.034d0 
ASTARN=1.0472d0 
BSTARN=0. 087266 
C 

ZETAN=  ( 2 .  0D0*ASTARN-  (RALN1+RALN2 ) )  /  (RALN2-RALN1 ) 

ETAN= ( 2 . ODO*BSTARN- (RBETN1+RBETN2 ) ) / (RBETN2-RBETN1 ) 

XN= ( 2 . 0D0*RAL- (RALN1+RALN2 ) ) / (RALN2-RALN1 ) 

YN=(2 . ODO*RBETA- (RBETN1+RBETN2 ) ) / (RBETN2-RBETN1 ) 

C 

ETI=  ( ( 5 .  ODO*  ( ZETftN**2 ) )  -  ( 4 .  ODO*ZETRN*XN ) -1 .  ODO  )  * 

+ ( ( (XN**2 ) -1 . ODO ) **2 ) / ( ( ( ZEIAN**2 ) -1 . ODO ) **3 ) 

C 

GN= ( ( 5 . ODO* ( ETAN**2 ) ) - ( 4 . 0D0*ETAN*YN ) -1 . ODO ) * 
+(((YN**2)-1.0D0)**2)/(((ETAN**2)-1.0D0)**3) 

C 

CNRB=M*FN*GN 

C 

IF  (RAL  .LT.  0.69813)  THEN 
C 

CNRB=0 . OdO 
GOTO  1000 
ENDIF 
C 

IF  ((RBETA  .LT.  -0.174532)  .OR.  (RBEIA  .GT.  0.34906))  THEN 
C 

CNRB=0 . OdO 
GOTO  1000 
ENDIF 
C 

1000  CMN= ( CMN1*EPA0  2  S ) + ( GNDAD*DDA )  +  ( ( CNDRD*DRUDD*DRFLX3 ) *EPA43 ) + 

+( (CNDTD*DTFLX3)*DELEDD)+(CMNP*PB)+(CMNR*RB)+(DCNB*BETA) 

++CNRB 

C 

(^**********  1HRUST  TERMS  ***********************************0 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 


n 


c 


c 


c 

c 


c 


c 


THIS  SECTION  DETERMINES  TOE  EFFECT  OF  THE  THRUST  VALUES  FOR 
ADDITION  TO  CX,  CY,  CZ,  OM,  CMM,  AND  CNM  VALUES  DETERMINED 
ABOVE  AND  CONTAIN  TOE  FOLLOWING  VARIABLES: 

CPTAL  -  COSINE  OF  PITCH  VECTOR  ANGLE 
SPTAL  -  SINE  OF  PITCH  VECTOR  ANGLE 
CYTAL  -  COSINE  OF  YAW  VECTOR  ANGLE 
SYTAL  -  SINE  OF  YAW  VECTOR  ANGLE 
ENGPQ  -  PORT  ENGINE  THRUST/ (QBAR*S) 

ENGSQ  -  STARBOARD  ENGINE  THRUST/ (QBAR*S) 

CXENGP  -  COEFFICIENT  OF  PORT  ENGINE  THRUST  IN  X  DIRECTION 

CXENGS  -  COEFFICIENT  OF  SBRD  ENGINE  THRUST  IN  X  DIRECTION 

CXT  -  COEFFICIENT  OF  TOTAL  THRUST  IN  X  DIRECTION 
CYENGP  -  COEFFICIENT  OF  PORT  ENGINE  THRUST  IN  Y  DIRECTION 

CYENGS  -  COEFFICIENT  OF  SRBD  ENGINE  THRUST  IN  Y  DIRECTION 

CYT  -  COEFFICIENT  OF  TOTAL  THRUST  IN  Y  DIRECTION 
CZENGP  -  COEFFICIENT  OF  PORT  ENGINE  THRUST  IN  Z  DIRECTION 
CZENGS  -  COEFFICIENT  OF  STARBOARD  ENGINE  THRUST  IN  Z  DIRECTION 
CZT  -  COEFFICIENT  OF  TOTAL  THRUST  IN  Z  DIRECTION 
CLMT  -  ROLL  MOMENT  COEFFICIENT  DUE  TO  THRUST 
CMMT  -  PITCH  MOMENT  COEFFICIENT  DUE  TO  THRUST 
CNMT  -  YAW  MOMENT  COEFFICIENT  DUE  TO  THRUST 

CPTAL=COS ( PTAL ) 

SPTAL=SIN(PTAL) 

CYTAL=COS ( YTAL ) 

SYTAL=SIN ( YTAL ) 

CRAL=OOS(RAL) 

SRAL=SIN(RAL) 

ENGPQ=ENGP/QBARS 
ENGSQ=ENGS/  QBARS 

CXENGP=ENGPQ*CPTAL*CYTAL 

CXENGS=ENGSQ*CPTAL*CYTAL 

CXT=CXENGP+CXENGS 

CYENGP=ENGPQ*CPTAL*  SYTAL 
CYENGS=ENGSQ*CPTAL*SYTAL 
CYT=CYENGP+CYENGS 

CZENGP=ENGPQ*SPTAL 
CZENGS = ENGSQ* SPTAL 
CZT=CZENGS+CZENGP 

CIMT=(CZENGS-CZENGP)*(25.5d0/12.0d0)/BWING 

OMT=CXT*  (0 . 25d0/12 .  OdO  )/CWING+ 

+  CZT*20 . 219d0/ CWING 

CNMT= ( CXENGP-CXENGS ) * ( 2  5 . 5d0/ 12 . OdO ) /BWING- 
+  CYT* 20 . 219dO/BWING 
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CX=CFZ*SRAL-CFX*CRAL+CXT 

Cy=CFY-t-CYT 

CZ=- ( CFZ*CRAL+CFX*SRAL ) +CZT 
CLM=CML+CIWr 
GMM=C«M+CMMT 
CNM=CMN+CNMT 
C 

C  THE  0.25/12.0  IS  THE  Z  OFFSET  OF  THE  THRUST  FROM  THE  OG 

C  THE  20.219  IS  THE  X  OFFSET  OF  THE  THRUST  FRCM  THE  OG 

C  THE  25.5/12.0  IS  THE  Y  OFFSET  OF  THE  THRUST  FRCM  THE  OG 

C 

C  RETURN  CX,  CY,  CZ,  CLM,  OM,  CUM  TO  CALLING  PROGRAM. 

C 

C 

RETURN 

END 
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